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SUMMARY 


The numerical calculation of unsteady two dimensional airloads 
which act upon thin airfoils in subsonic ventilated wind tunnels is 
studied. Neglecting certain quadrature errors, Bland's collocation 
method is rigorously proved to converge to the mathematically exact 
solution of Bland's integral equation, and a new three-way equivalence 
is established between collocation, Galerkin's method and least squares 
whenever the collocation points are chosen to be the nodes of the 
quadrature rule used for Galerkin's method. The computer program 
displays remarkable convergence with respect to the number of pressure 
basis functions employed, amd agreement with known special cases is 
demonstrated. New results are obtained for the combined effects of 
wind tunnel wall ventilation and wind tunnel depth to airfoil chord 
ratio, including acoustic resonance between the airfoil and wind 
tunnel walls. 
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but it will be clear from the context which is correct. Any consistent 
system of physical units may be employed, although most quantities are 
nondimensional . 
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jl- Introduction 1 


Recent interest in aerodynamic testing at high subsonic and tran- 
sonic speeds has increased the need for improved understanding of the 
phenomena involved, and for improved computational methods to guide 
experiments as well as to serve as a rational means for extrapolating 
wind tunnel test data to free flight conditions. 

While fundamental differences exist between two and three dimen- 
sional flows, solutions of two dimensional problems provide advance 
insight into three dimensional flow phenomena. Calculations for two 
dimensional problems are significantly simpler, and testing can be 
performed using models with constant sections extending completely 
across the tunnel, or enclosed between splitter plates. 

Aerodynamic interaction between the wind tunnel walls and the model, 
commonly termed interference, is present in all subsonic wind tunnels 
and becomes more complicated at higher subsonic and transonic speeds. 
Efforts to utilize opposing interference effects associated with closed 
wall and open jet tunnels have resulted in the development of ventilated 
wind tunnels at various facilities. Whereas theories of interference 
are highly developed for steady subsonic flow in fully closed or open 
jet tunnels, relatively little is known about unsteady flow, especially 

■ • 9 Q 

in ventilated tunnels. 


1( rhe authors wish to acknowledge the help of Dr. Sanford S. Davis, 

Ms. Sara Alfont, Messrs. Tuli Haromy, Charles Doughty and Karl Kuopus. 

2 A comprehensive survey of the state of the art may be found in the work 
by Garner, Rogers, Acum and Maskell [1,1966], See also Goethert [2,1961], 
Glauert [3,1933] and Pope and Harper [4,1966]. 

3 Numbers in square brackets refer to the bibliography in the order of 
citation, followed by the year of publication, and possibly the page, 
section or chapter of interest. 
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This work addresses the problem of predicting unsteady airloads on 
oscillating thin planar airfoils in subsonic ventilated wind tunnels, and 
extends the earlier work of Bland [5,1970] which is based on an approximate 
ventilation boundary condition that is exact for only the closed wall and 
open jet conditions. 

We rigorously establish in Section 4 that Bland's collocation method 
must converge to the mathematically exact solution of his integral equation. 
The analysis is made using two sets of orthogonal polynomials, one for 
downwash the other for pressure. These polynomials, called airfoil poly- 
nomials, enable an elegant formulation in terms of dual generalized Fourier 
series. The known closed form Sohngen and Kussner-Schwarz solutions are 
reformulated in terms of airfoil polynomials, and used to evaluate the 
accuracy of the computer calculations. 

The computer program developed during this study is called TWODI and 
will calculate pressures, section coefficients and generalized aerodynamic 
forces for arbitrary combinations of Mach number, reduced frequency, tunnel 
depth to airfoil chord ratio, ventilation coefficient and downwash. TWODI 
is shown to be in accurate agreement with known closed form results and 
extends previous results of unsteady flow to ventilated tunnels. Various 
calculations are presented, including predictions of the combined effects of 
depth to chord ratio, and of acoustic resonance between the wind tunnel walls 
and the airfoil. 

While the reader who wishes to avoid theoretical considerations cannot 
help but be at some disadvantage, it is possible to use the TWODI program 
by proceeding directly to Section 7. Other sections may then be read as 
the interest arises. 


§2. History of the two dimensional problems 


This section presents an extensive although not complete history 
of developments in the problem of predicting airloads in unsteady two 
dimensional subsonic flow. Attention is given to the emergence of 
exact closed form solutions because these provide permanent and en- 
during standards of comparison for numerical as well as theoretical 
considerations. The integral equations relating downwash and pressure 
in linearized potential flow are especially emphasized. 

Three dimensional results are generally excluded except insofar as 
they relate to the two dimensional problem. Nevertheless we note that 
one of the most important uses of two dimensional theories is to provide 
insight into the practical matters of three dimensional problems. 

The most important engineering application of unsteady aerodynamics 
traditionally has been to calculate forcing functions for structural 
dynamics and aeroelasticity. This interest is readily apparent in the 
earliest papers, and it may be observed that a recurring practical 
problem has been to maintain acceptable accuracy at higher frequencies. 

While the physical conditions which justify linearized inviscid theories 
are not satisfied at arbitrarily high frequencies, an arbitrary frequency 
capability is worthy because the ability to analyze arbitrary time depen- 
dence demands it mathematically. 

Vfe begin with a brief description of the airfoil equation. By 1918, 
Prandtl at Gottingen had completed the theory of bound vortices, and 
Ackermann had performed calculations for steady lift on airfoils at low 
speed. However, Ackermann* s work was interrupted by the war and was not 
reported until afterward by Birnbaum [6,1923], then only at Prandtl* s 
request. Comparable calculations were performed in the United States 
by Munk [7,1922], [8,1924]. Restated in terms of pressure jump rather 
than vorticity, it was necessary to solve an integral equation of the form 

1 

w(x) = / K(x~£)Ap(£)d£ (2- 

where the kernel is given by 
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and subject to the Kutta condition that the pressure jump vanish at the 
trailing edge 

lim Ap(£) = 0 1 . (2-3) 

The downwash, assumed known, is given for steady flow by the streamwise 

derivative of the vertical coordinate of the airfoil contour: 

, . dh .. 

w(x) = — • (2-4) 

The coordinate system and airfoil contour are shown in Figure 1. 



Figure 1. Coordinate system and airfoil contour 

We shall refer to the integral equation (2-1) with kernel given by 
equation (2-2) and subject to the Kutta condition (2-3) as the airfoil 
equation . It may be observed that the airfoil equation is a singular 
Fredholm integral equation of the first kind and its kernel is of difference 
type with a Cauchy singularity. 2 Much of the difficulty encountered in 


^■We shall generally, but not always, follow the practice of stating values 
of the downwash function in terms of x and values of the pressure function 
in terms of £ for the purpose of emphasizing certain symmetries. 

2 An integral equation of the form 

f(x) = f b K(x,?)u(£)d£ 
cL 

where f is a known function and K is a known function (called the kernel 
function) is said to be a Fredholm integral equation for u of the first 
kind. An integral equation of the form 

f(x) = U(x)+/ b K(X,?)u(C)d? 

3 . 

is said to be a Fredholm integral equation of the second kind. If the 
kernel depends upon the difference x-5, it is said to be of difference 
type and we write K(x-£) . If the kernel is bounded for all values of x 
and the integral equation is said to be nonsingular, otherwise singular. 

If the kernel is of difference type and of the form 

K(x-?) = x^T 

where C is a constant, it is singular and is said to have a Cauchy singularity. 



c o 


Tt4 ~JU 
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solving the airfoil equation and related equations in aerodynamics may 
be attributed to the fact that until recently, no general theory existed 
for the solution to integral equations of the first kind. The transla- 
tion and publication of Ivanov's book [9 , 1976] indicates that substantial 
activity has been taking place in the Soviet Union in this area, much of 
it unreported in the Western literature. Without detailed examination 
of this work its relationship (if any) to ours remains unclear - 

The airfoil equation, based upon classical linear potential flow 
theory, is rigorously valid for zero Mach number steady two dimensional 
flow about an airfoil with profile differing infinitesimally from the 
x-axis. Noting the linearity of the airfoil equation, Birnbaum calculated 
the first order effect of thickness by using the difference between the 
upper and lower surface profiles, and calculated the lift and moment using 
the centerline profile, for profiles described by third degree polynomials. 
The pressure functions 



when substituted into the airfoil equation produce downwash functions that 
are constant, linear and quadratic, respectively, and satisfy the Kutta 
condition. The Birnbaum-Ackermann calculations consist of integrating 
in closed form each of the above basis functions (Grundfunktionen) to 
obtain an induced downwash (induzierte Vertikalgeschwindigkeit) , which 
is integrated in turn to produce a profile function. These profile functions 
are then superimposed to represent an arbitrary profile of third degree. 

The general procedure is to represent the pressure as a linear combination 
of NP prescribed pressure basis functions 

Ap(£) = a : /^|+ a 2 /l-£ z + ... + ^l-£ 2 (2-5) 

and to determine the unknown coefficients according to a prescribed profile 
which is given by a corresponding linear combination of basis functions, 
in this case, powers of x, 

h(x) = b 2 x + b 2 x 2 + ... + bjjpX 1 ®. (2-6) 

Although no discussion was made of collocation, residual error or con- 
vergence, the method of solution appearing in the Birnbaum-Ackermann paper 
of 1923 is a forerunner of modern collocation solutions to the integral 
equations of aerodynamics. 
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The first published calculations for unsteady flow are due to 
Birnbaum [10,1922] and are closely connected with the Birnbaum-Ackermann 
paper. Birnbaum was a doctoral student of Prandtl's and his dissertation 
was devoted to extending Prandtl's theory of bound vortices to the more 
difficult problem of unsteady flow. Through an analysis of free and 
bound vortex sheets, Birnbaum was led to infinite series of the form 

E k n ( log k) n (2-7) 

m>n 

which did not converge well even for small values of reduced frequency 
near k = 0.10. Birnbaum apparently coined the term reduced frequency , 3 
and defined it in terms of the semichord as 

k = — • (2-8) 

v 

CO 

He noted that the numerical limitation 

k < .1 

corresponded physically to less than one structural oscillation per thirty 
chordlengths of flow, and referred to this as a quasi-steady oscillation 
(quasi-stationare Schwingungen) - Birnbaum analyzed time dependence with 
complex Fourier series and introduced into unsteady aerodynamics the con- 
cept of complex amplitude of motion 


h(x,t) = Re[h(x)e lut ] 


(2-9) 


for the particular cases of plunging and pitching, but did not use complex 
downwash nor complex pressure 


w (x , t ) = Re [w(x) e itijt ] = Re[(^- + ik) h (x) e itot ] , 


( 2 - 10 ) 


Ap«,t) = Re[Ap(?)e iut ]. (2-11) 

Lift and moment were graphed as complex vectors parameterized by reduced 
frequency from k = .00 to k = .12. Birnbaum’ s interest in structural 
dynamics was clear. He treated the static and dynamic stability of an 
elastically suspended wind tunnel model using the concept of energy trans- 
fer and presented damped and divergent experimental response curves 
together with theoretically calculated frequencies. 

3,, reduzierte Frequenz", Birnbaum [10,1922, p,12] , [11,1924, p.279]. 
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It is significant that Birnbaum did not find the integral equation 
between pressure and downwash for unsteady incompressible flow. This 
equation is the same as the airfoil equation (2-1) except that the down- 
wash and pressure are complex as in (2-9) and (2-10) , and the kernel 4 
is given by 

K(x,k) = -j— - e“ ikx [Ci(k|x| ) + iSi(kx) + (2-12) 

where Ci and Si are the cosine and sine integrals. 5 This equation was 
not derived until 1938 by Possio as a limiting case of what we now call 
Possio' s integral equation for compressible subsonic unsteady flow. 

The next work on unsteady flow after Birnbaum 1 s was published by 
Wagner [20,1925 1 6 based on his doctoral dissertation at Berlin under Hoff 
and Hamel. Wagner used the method of conformal mapping to solve the two 
transient problems of a flat plate accelerating from rest to constant 
velocity, and of a step change in angle of attack. Glauert [21,1929] used 
the same method to compute lift and moment on an oscillating airfoil for 
reduced frequencies up to k = .5. Kussner [22,1929] extended the frequency 
range further to k = 1.5, avoiding the slowly convergent series 7 of Birnbaum 
and expressing the pressure as 

Ap(?) = + a l A-C 2 + a 2 C A-C 2 + a 3 (4£ 2 -l) A-? 2 + ... (2-13) 

By making the coordinate transformation, 

e = -cos 0, (2-14) 


Kussner was able to simultaneously recast the pressure basis functions 
as a Fourier sine series except for the first term 

0 °° 

Ap(0) = agcot-j + a n sin n0 (2-15) 


4 Seev e.g. , Bland [12,1968,54.1], [5,1970,18]. 

5 Se^ e.g., Abramowitz and Stegun [13,1964,15.2]. The complex setting in 
Abramowitz and Stegun does not render Ci an even function. To avoid error 
in the logarithmic term, the absolute value sign is necessary in (2-12) , 
and in this regard, it could be considered as missing or at least ambiguous 
in Possio [14 , 1938 ,p. 448] , Dietze [15 , 1946 ,p. 22] Watkins, Runyan and 
Woolston [16,1955, App. B] Bisplinghof f , Ashley and Halfman [17 , 1955 ,p. 325] , 
Garrick [18 , 1957 ,p . 710] , and Fung [19, 1969, p. 429] . 

6 For further discussion, see, e.g., Y.C. Fung [19, 1969 , § §5 . 8 , 6 . 7 , 15 . 1] or 
Bisplinghoff , Ashley and Halfman [17 , 1955 , § 5 . 7] . 

7 See, H.G. Kussner [23,1953, §3] for further information. 
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and the downwash as a Fourier cosine series 

w ( 0 ) = bn + 2 ? b n cos n0 . (2-16) 

n= i 11 

He analyzed the oscillating flat plate and approximated the discontinuous 
downwash of a flap with a fifth degree polynomial. 

Theodorsen [24,1935] obtained the first closed form solution for the 
force and moment on a flat plate undergoing pitch rotation, vertical trans- 
lation and aileron rotation. Theodorsen 's solution, though restricted to 
rigid body motions, was valid for an arbitrary frequency range. He used 
the method of conformal mapping employed earlier by Wagner. The flutter 
solution for the three degree of freedom system was thus reduced to a non- 
Hermitian three by three matrix eigenvalue problem, and the three special 
cases called torsion-aileron, aileron-deflection, and def lection-torsion 
were solved in closed form, producing flutter velocity curves vs. in vacuo 
frequency ratios. Theodorsen improperly expressed his intermediate results 
in terms of a quotient of nonconvergent integrals 


c(k) 


a -ik£_ 


/s 2 -i 


d? 


oo 




e 


-ik£ 


_£+l_ 
/? 2 + 1 


which he identified in terms of Hankel functions of the second kind 


C (k) = 


Hi ( 2) (k) 


Hp } (k) + 1 h{ 2) (k) 


(2-17) 


where 

Hp } (k) = J n (k) - iY n (k) (2-18) 

and where Jn and Y n are standard Bessel functions of the first and second 
kinds. 8 The function C(k) is called Theodorsen 1 s circulation function. 

Although his derivation of it was faulty, the result was shown later by 
Schwarz to be correct. 8 


8 See, e.g., Abramowitz and Stegun [13,1964]. 

8 Schwarz [25,1940,14]. see also the discussion in Bisplinghof f , Ashley 
and Halfman [17 , 1955 ,p. 272] . 
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Shortly after Theodor sen's solution in 1934 for three downwash functions, 
the same problem for arbitrary downwash functions was solved independently 
by Cicala [26,1935] and Kussner [27,1936], They showed that the coefficients 
of the downwash and pressure basis functions appearing in equations (2-15) 
and (2-16) are related recursively: 


a 0 = 


_ b 0 H| 2 > (k) + ibiH^ 2 ) (k) 

h{ 2) (k) + i H^ 2 ) (k) 


ik , , ik , ^ _ 

a n = ^ b n-l “ bn “ 2n bn+i; n - 1 ’ 


(2-19) 


(2-20) 


The importance of the Cicala-Kussner extension to arbitrary downwash 
functions is notable. While rigid body motions as discussed by Theodorsen 
are important, higher flexible modes are essential to the theory of aero- 
elasticity, Whereas chordwise flexibility is typically insignificant in 
two dimensional flow, it has long been recognized as significant in the 
design of modern airplanes , 1 0 and the ability of a two dimensional aerodynamic 
method to shed practical light on three dimensional methods requires that it 
possess an arbitrary chordwise downwash capability. 

The problem of unsteady airloads on an airfoil entering a sharp-edged 
gust was solved for incompressible flow by von Kdrmdn and Sears [29,1938] 
using conformal mapping. In the same year, Possio [14,1938] obtained the 
kernel relating the downwash and pressure for compressible subsonic flow 
about an oscillating airfoil: 


K(x, k,M) = - e" lkx {exp(^^) [iM sgn(x)Hp ) (^z^-) + H(j) 2) (^j 

L + ik /* exp(^r)Hj 2 > (Js^jjr^dA}. 


+ — e log 

TT 


( 2 - 21 ) 


Although Possio did not explicitly write down the integral equation, 11 he 
performed numerical calculations for lift and moment on an oscillating flat 
plate at values of reduced frequency up to k = .6 and at four values of 
Mach number, M = 0,.25,.50, and .70. Using the same pressure basis 


10 See, e-g.. Turner, Clough, Martin and Topp 


[28, 1956, p.805] . 


11 The Possio integral equation first appears in Kussner' s Allgemeine 
Tragflachentheorie [30,1940, §7] as a special case of the Kussner integral 
equation for three dimensional flow. 
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functions as Birnbaum, Possio observed that the solution to the integral 
equation required that the linear combination of downwashes induced by 
the pressure basis functions must equal the prescribed downwash at all 
points on the airfoil? i.e., 

NP 

w(x) = I a n w n (x) , (2-22) 

n=i 

which is in general not possible- Possio resolved this difficulty by 
collocating (2-22) at an equal number of discrete points according to the 
system of linear algebraic equation 


NP 

n | w n < x m > a n = w(Xm) ? m = 1, . . . ,NP. (2-23) 

This is the earliest use we have found of the method of collocation to 
solve an integral arising in aerodynamics. 

We note that a least square error solution of equation (2-22) in which 
the overall error is integrated numerically and then minimized results in 
the matrix equation 

[w n (Xm)]*[ W m ] [w n (x m )]{a n } = [w n (x m )]*[ W m ]{w(x m )}, (2-24) 

NPxNQ NQxNQ NQxNP NPxl NPxNQ NQxNQ NQxl 


where x m ; m = 1,...,NQ are the nodes of the quadrature rule and 
W m ; m = 1 , . . . , NQ are its weights, and where * denotes the complex con- 
jugate of the transposed matrix, a fact sometimes missed. 12 

Under fairly restrictive assumptions, Sohngen [32,1939] proved that 
the solution to the airfoil equation is given by 


Ap(£) 


4 / l-g 1 / l+x 

7T / 1+5 L X / 1-X 


w(x) 

C-x 


dx. 


(2-25) 


which is called the Sohngen inversion formula and is a particular case of 
a more general formula given in [33, 1953, Ch 11]. These restrictions were 

later weakened by Sohngen [34,1954] to include all downwash functions of 

, ~ 4 

class LP,p>l. J A weaker proof based on the stronger assumption that p>-j 

was given by Tricomi [35,1951] and appears in his book on integral equations 


12 Since the error is a real valued noncoristant function of several complex 
variables, it is not differentiable. A rigorous proof of (2-24) is given 
by Fromme [31,1964, Appendix A], 

13 A function f belongs to class LP over the interval I if /Jf (x) |^dx is 
finite. 
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[36,1957], Sohngen' s later proof and Tricomi's proof both identify the 
integral equation (2-1) with the difference kernel (2-2) as a finite 
Hilbert transform. 14 We point out that these mathematical results make 
essential use of the condition that the pressure be finite at the trailing 
edge (and hence zero, thereby satisfying the Kutta condition (2-3)), for 
otherwise the finite Hilbert transform has no unique solution in the iP 
space, p>l, and contains an additional term of the form 

,C . 


where C is an arbitrary constant. Therefore, if the Kutta condition is 
denied, the solution to the airfoil equation is not unique. 

tt 

Five important papers appeared in 1940: (1) Kussner published the 

integral equation relating pressure and downwash in unsteady three dim- 
ensional flow; 15 (2) Kussner and Schwarz [39] , [40] analytically summed 
the series (2-15) , (2-16) , (2-19) and (2-20) to provide in closed form 

the pressure in terms of downwash for oscillatory flow as 


A 4 r i r r sm 0 . lk . , 1-cos (n+6) 

Ap (6 ) = — f { + — sm n log — 

tt 0 cos 0-cos n 2 l-cos(n-6) 


+ cot — [l+cos n+d-cos n)T(k) ] }w(n)dn 


(2-26) 


where 


T(k) 


Hp> (k) - i Hp> (k) 

Hp> (k) + i Hp' (k) 


= 2C(k) 


1; 


(2-27) 


II 

(3) Kussner [41,1940] used the method of Laplace transforms to extend the 
solution (2-25) to arbitrary time dependence? (4) Schwarz [25,1940] used 
the Sohngen inversion formula to obtain the solution (2-26) by an independent 
approach; and (5) Sohngen [42,1940] used the inversion formula (2-25), 


14 The finite Hilbert transform of cf) is the function f with values given by 


f (x) 


1 r l <t> (£) 

IT ”1 £-x 


d£. 


See also Sneddon [37 , 1971 ,p. 238, p. 467] . 

15 H.G. Kussner [ 30, 1940], [38 , 1940] . However, the kernel was not put into 
a computationally tractable form until fifteen years later by Watkins, 
Runyan and Woolston [16,1955]. 


f 
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together with Birnbaum's method and Laplace transforms to derive the 
solution for the oscillating airfoil and for an airfoil starting from 
rest. Equation (2-26) is called the Kussner-Schwarz solution and in 
Cartesian coordinates it is given by 

Ap(U = £ + l-C(k)]-ikA(x,£) /jzj* w(x)dx, (2-28) 


where 


A(x,£) = log 


/i±i+ /Al* 

/ 1-g / 1— x 


m-/¥- 


+x 


(2-29) 


Taking advantage of the theoretical developments principally by 
Theordorsen and Kussner, Smilg and Wasserman [43,1942] published their 
AFTR 4798 which was to serve for years as the standard reference in the 
United States on flutter. The kernel to Possio's equation was tabulated 
by Schwarz [44,1943], and Possio's numerical calculations based on the 
method of collocation were extended by Fraser [45,1941], Fraser and Skan 
[46, 1942] , and Schade. 16 Dietze, 17 working under the supervision of 

ii 

Sohngen, introduced an iteration method to solve Possio's integral equation 

ii 

approximately, making use of the known Kussner-Schwarz solution. Dietze' s 
iteration method consists of, given the downwash, calculating the pressure 
according to the Kussner-Schwarz solution for incompressible flow, and 
calculating the downwash required by that pressure according to Possio's 
integral equation. The residual error in downwash furnishes the starting 
point for successive iterations. Dietze* s calculations were for values 
of reduced frequency up to k = 1.0 for pitch rotation, vertical translation 
and flap rotation. In a major review of the two dimensional problem, done 
in two parts by Karp, Shu and Weil [52,1947], and Karp and Weil [53,1948], 
it is observed that in the case of a discontinuity in downwash the Dietze 
iteration will not eliminate the logarithmic singularity at the discontinuity 
so that slow pointwise convergence might be found in the neighborhood of 
a discontinuity. Additional calculations with Dietze' s method by Turner 
and Rabinowitz [54,1950] indicated that numerical values of overall section 


16 T. Schade [47,1944]. An English translation is available as [48,1946]. 

17 F. Dietze [49,1943], [50,1944], English translations are available as 

[15,1946] and [51,1947]. A description of Dietze 's method may be found in 
Y.C . Fung [19,1969, §14. 5] . 
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coefficients appeared to be acceptable, but that the rate of convergence 
deteriorated for larger values of Mach number and reduced frequency. 

Apparently no mathematical proof of convergence of Dietze 1 s method has 
ever been given. A compilation of tables from various references was 
given by Luke [55,1950], and another method which approximated the continu- 
ous part of Possio's kernel by a polynomial was introduced by Fettis [56,1952], 
The next major advance comparable to the Kussner-Schwarz solution is 
a general solution to Possio's integral equation using elliptic coordinates 
and Mathieu functions. It is obtained by starting from the basic differ- 
ential field equations and boundary conditions 

v 2 <j) - M 2 (— — I- ik) 2 <)> = 0 (2-30) 

OX 

P = -2(^- + ik)*, (2-31) 

w(x) = lim 1^-, |x| < 1, (2-32) 

y+0 9y 

H IP 

Upon making the Kussner transformation in the form 10 


- - y - k kM 

x = x, y = -g, k = -gj-, k = -g?, 


(2-33) 


<j > = Be-*- 10 ^, w = Be ; ‘- I<x w, p = e iKX p, 


(2-34) 


there results 


V 2 <j> + k 2 4> = 0, (2-35) 

P = -2(^4 + ik <f) , (2-36) 

dX 

w (x ) = lim -^4, I x I <_ 1. (2-37) 

y->o <>y 

Applying Green's identities, integrating (2-36), and introducing elliptic 
coordinates (£,n) according to 

x = -cosh f cos n, y = sinh £ sin n# (2-38) 


18 This transformation, valid for arbitrary time dependence in three dimen- 
sional flow, was first used in aerodynamics by Kussner [30,1940, §3] , and is 
a combination of the Galilean and Lorentz transformations. See, e.g.. Miles 
[57,1959, §§ 2. 3,2.4] . For the special case of steady flow, Kussner* s trans- 
formation reduces to the Prandtl-Glauert transformation. 
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Kussner [ 2 3, 1953, § §2-4] showed that 

p(£,6) = -|- /g { G,n(K,?,o,0) [g, 0 (K,o,n,Tr)T(K,k)-G, e (K,o,n,o) ] 

+ (ik — r~ — ^-)G(K,C,n»0) }w(n) sin n do, (2-39) 

sm n on 

where 

Ne^ 2) (^-,C) 2 2 

G (<,5,0,0) = £ \ se n (^~, 0 ) s e n (^r~, 0 ) , (2-40) 

11=1 Net2)(j£^7 4 4 

and 

co-in / 2 s\ 

£co+i w /2 ex P (_ik cosh ?)G, n e(K,C,0,ir)d? 

T(K ' k) = • < 2 - 41 > 

-oo+iir/2 eXp(-ik COsh S) G 'n0 (K*?»°/0)dC 

In the above, partial differentiation is denoted with a comma followed by 
a subscript, and se n and Ne^ 2 ^ denote Mathieu functions of the first and 
third kind, respectively. 19 Equations (2-39), (2-40) and (2-41) reduce in 
the case M = 0 to equations (2-26) and (2-27) . 

The application of elliptic coordinates and Mathieu functions to the 
unsteady flow equations was undertaken first by Reissner and Sherman [59,1944] , 
followed with apparently independent work by Biot [60,1946], Timman [61,1946], 
and Haskind [62,1947]. 20 Among these, Timman ’s was the most successful. 

Further work was done by Billington [64,1949], Timman and van de Vooren 
[65,1949], Reissner [66,1951] and [67,1951], Timman, van de Vooren and Griedanus 
[68,1951], Kussner [23 f 1953], 21 van Spiegel and van de Vooren [69,1953], 
de Jager [70,1954], and Williams [71,1955]. Despite its success, the method 
of elliptic coordinates and Mathieu functions has been slow in gaining the 
recognition it merits. 

Using successive conformal mappings and Jacobian elliptic functions, 

Timman [72,1951] obtained in closed form the exact solution for pressure on 


19 See, e.g., N.W. McLachlan [58,1947]. 

20 In Russian. An English translation is available as [63,1947]. 

21 This contains a comprehensive and elegant discussion of the two dimen- 
sional problem up to that time. 
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an oscillating airfoil in a wind tunnel with parallel closed walls. 

Tinman' s solution is for incompressible flow and reduces to the Kussner- 
Schwarz solution in the limiting case of infinite tunnel depth. The analysis 
applies to a centrally located airfoil as shown in Figure 2. 



Figure 2. Airfoil in a wind tunnel 

The condition that the tunnel wall be closed is the kinematical one that 
the component of velocity normal to the wall vanish at the wall. 

jr <f> (x,±n H ) = 0. (2-42) 

For incompressible flow, wind tunnel wall effects tend to become 

maximum at some small value of reduced frequency and to diminish as the 

frequency is increased further. Disturbances produced by an oscillating 

airfoil are reflected by the walls instantaneously in the sense that the 

acoustic transit time for a given distance is much less than the flow transit 

time for the same distance. Mathematically this ratio of transit times is 
v 

zero because — = M = 0. However, when the fluid is compressible, a nonzero 
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transit time is required and it is possible for the frequency of airfoil 
oscillation to be such that succeeding disturbances arrive back at the 
airfoil so as to reinforce one another, thereby causing an acoustic 
resonance. Runyan and Watkins [73,1953] predicted that this resonance 
phenomenon would occur at frequencies given by 

w n = 7 ~(n-l); n = 1,2,... (2-43) 

These resonant frequencies are finite for all subsonic Mach numbers in a 
compressible fluid, and become lower as the Mach increases. The resonant 
reduced frequencies are 


kn " Mrljj (n 2 )? n ' 1 ' 2 '"- 


(2-44) 


and decrease rapidly at higher Mach numbers. For example, 

/3 

M = — & n H = 10 -* k x = .090690. 


(2-45) 


Thus, acoustic resonance between the airfoil and the tunnel becomes an 
important phenomenon at high subsonic speeds because it occurs at relatively 
low values of reduced frequency. Furthermore, these resonant frequencies 
are the same for all down washes. 

Runyan and Watkins also showed that the pressure and downwash are 
related for subsonic compressible flow in a closed wind tunnel by an 
integral equation of the form (2-1) where the kernel may be expressed as 
the sum of the Possio kernel (2-21) plus an incremental part representing 
the presence of the walls. The incremental kernel becomes singular at the 
resonant frequencies (2-44) and since the kernel may be interpreted physically 
as the downwash at one point due to unit pressure at another point, the 
pressure should drop to zero at resonance. This effect is entirely analogous 
to the elementary spring-mass system in which vanishingly small forces can 
produce large amplitudes when the excitation frequency is sufficiently 
close to the resonant frequency. 

The above theoretical predictions were verified experimentally and 
computationally by Runyan, Woolston and Rainey [74,1956] . The kernel 
function was approximated by assuming the tunnel depth to chord ratio to 
be large, the pressure was represented by the first three terms of the 
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series (2-15) , and the integral equation was collocated at three points. 

Good agreement was obtained between theory and experiment for phase angles 
and fair agreement was obtained for magnitudes. Resonant frequencies were 
predicted accurately by the theory, and the expected drop in lift at 
resonance appeared to be quite abrupt. 

Bland, Rhyne and Pierce [75,1967] extended the theory to include narrow 
channel flow in connection with destructive oscillations in nuclear rocket 
engines. The kernel obtained previously by Runyan and Watkins was expressed 
as a single function in a simpler form without using the deep tunnel 
approximation : 


(x,k,M,n H ) = ~*(l+sgn x)e ikx tanh knn 


r 7rR n x , 

ikM 2 x f sgn x ik GXp BriH 

+ “’“’‘HF-’ni, ^ + ,-S-bL,? 


where 


*n = /(n-i ) 2 - <^> 2 . 


(2-46) 


(2-47) 


At the resonant frequencies given by (2-44) , the kernel (2-46) becomes 
infinite because R n = 0. The kernel is singular at x = 0 and the infinite 
series in (2-46) behaves in the neighborhood of 0 like a slowly convergent 
geometric series. The pressure was represented in the form 

sin[(n-y)cos _1 £] 

1 ( 

but no elaboration was offered for the reason behind this choice. Never- 
theless, for continuous downwash functions, this marks the first basic 
improvement in the choice of pressure basis functions since 1929, 22 and 
as we shall see, leads to an elegant closed form generalized Fourier 
solution for pressure in certain special cases. 

Because the physical problem of interest to Bland, Rhyne and Pierce 
involved low speed flow, numerical calculations were performed only with 

22 Cf . equation (2-15). The singularities in pressure that must accompany 
downwash discontinuities have been given by White and Landahl [76,1968] , 
Land ah 1 [77,1968] and Rowe, Sebastian and Redman [78,1976]. 


/ 1 -r oo 
= / l+E n= 
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M = 0 using a form of the kernel (2-46) in which certain of the series 
could be summed in closed form. The integral equation was collocated 
using Hsu's technique [79,1958] in which collocation points are inter- 
digitated with the quadrature points so as to minimize the error in lift 
on the average . Flutter calculations using two degrees of freedom were 
performed for depth to chord ratios from 0.06 to 2.00 and verified by 
experiment. 

The problem of unsteady flow in ventilated wind tunnels was under- 
taken by Bland in his doctoral thesis [12,1968] and subsequently published 
in the open literature [5,1970]. Bland approximated the boundary condition 
at a ventilated wall by a linear combination of the boundary conditions for 
open and closed walls. 


P + 


c w If = 0 at 


(2-49) 


where the semi- empirical constant c w is called the wall ventilation co- 
efficient (D. Davis and D. Moore [80,1953] discuss an approximation of c w 
for a wall with longitudinal slots.). The boundary condition (2-49) is 
exact only for the two limiting cases; (1) c w = 0, which corresponds to 
an open jet and, (2) c^ = 00 , which corresponds to a closed wall. Treating 
the pressure as an odd function of y for a centrally located airfoil. Bland 
applied the method of Fourier transforms to derive the kernel for the 
integral equation between downwash and pressure. Bland's kernel can be 
expressed as 


K(x,k,M,n H ,c w ) = _ -AjL log | x | + 


1+sgn (x ) 1+cyjk tanh kriH e -ikx 
c W + i” tanh kriH 


1 r / % _ a >|x|. ikn H IxL , ikM 2 x 

- T— [sgn(x)F' (^— *-) - F (*r-~) ] exp ^2 

prifj p p^h p 

* ^ - H? * <-* ■“» 3S5? 

ik [log tanh -757-) + (exp - -1) log tanh 1 , 


4tt£ 


4 Bn* 


43hh 


(2-50) 
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where 


F(6) = E {^2- exp(-X n 6) exp[-(n-i)ir6]}, 

n= l Xn A 2 


77 (n-j) 


(2-51) 


[1 + 


l+(YX n )‘ 


A n 

(2-52) 

MkriH 

(2-53) 

ex n ' 

0, 

(2-54) 


(2-55) 


t3.ii Xj^ + yX n — 0/ 
c w 

Y tih’ 

Equations (2-50) to (2-55) are slightly modified from the form originally 
given. Bland used two functions in place of F and F' , and we have identi- 
fied one of them as the derivative of the other, resulting in a numerical 
simplification. The series (2-51) has been accelerated by subtracting 
asymptotic terms which can be summed analytically. Even with this im- 
provement, due to Bland, the series are poorly convergent for very small 
values of the argument, especially for unsteady flow. 

The singularities of the kernel are confined to the first two terms. 
The first term, a Cauchy singularity, is the dominant one and is present 
in all two dimensional subsonic kernels. The second term is a weaker, 
integrable singularity of logarithmic type and is present in all the 
unsteady kernels. The remaining terms are bounded and can be accurately 
integrated using an appropriate Gaussian quadrature rule. Bland's complete 
kernel reduces in special cases to all kernels given previously. 

The method of solution employed by Bland is unique and is particularly 
interesting. From the Sohngen inversion formula. Bland observed that 
the function 


p « ) = / A § Ap(c) 


(2-56) 


results in an integral equation of the form 


w(x) = / 
-1 


1-C 

i+e 


K(x-£ ) P (£ ) d£ 


(2-57) 
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and that in the special case of the airfoil equation if w is a polynomial/ 
then P is a polynomial of the same degree . This property is based on the 
Bland integral transforms and will be discussed further in Section 3. Using 
this result as a starting point for equations with more complicated kernels , 
Bland introduced airfoil polynomials and used them to evaluate integrals 
of the singular part of the kernel in closed form and integrals of the 
bounded part using Jacobi-Gaussian quadrature. Collocation and quadrature 
points were interdigitated following Hsu's technique, and remarkably good 
numerical convergence with respect to the number of pressure basis functions 
was obtained. Numerical results including flutter calculations were presented 
for several problems of interest. 

Of the methods considered in this study, Bland's is the most germane 
to the problem of unsteady two dimensional subsonic flow in ventilated 
wind tunnels, it has superior numerical characteristics, and it is the 
most general. 

Convergence of Bland's method is proved for the first time in Section 4. 
We note that the proof of convergence of Bland ' s method hinges on the Cauchy 
singularity and is largely independent of the particular form of the remaining 
part of the kernel. 
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§3. Theory of airfoil polynomials 1 

The theory of airfoil polynomials springs from the solution to the 
airfoil equation for steady incompressible flow using the Bland trans- 
form, and may be applied to the numerical solution of any problem for 
which the linear transformation from pressure to downwash is a singular 
integral equation of the first kind whose kernel is dominated by a 
Cauchy singularity. 

Let 

cos [ (n-y) cos _1 x] 

X n < x > = y 3 ; n = 1,2,... (3-1) 

COS[yCOS 1 X] 
sin [ (n-y) cos” 1 £ ] 

<Pn(Z) = 1 ? n = 1 , 2 ,... (3-2) 

sin [ycos ~ 1 £ ] 

be defined on the open interval (-1,1). These functions, shown in Figure 3, 
are polynomials of degree n- 1 : 

Xi (x) = 1, X2< x > = “l+ 2 x,..., x n +2 (x) = 2 xx n +i ( x )“X n ( x ) ' ••• ( 3 “ 3 ) 

ih(t) = 1, 4*2 ) = 1 + 2 ^,..., ip n+ 2 (Z) = 2 £i |/ n+1 (C)-* n (C) f ... ( 3 -4) 

The above recursion formulas are neutrally stable so that computational 
errors are neither damped nor magnified with increasing n. 

These polynomials are orthogonal with respect to reciprocal weight 
functions 

1 1 ' I -|_v 

- Xm(x)x n (x)dx = 6 ^, (3-5) 

7 = 6^. (3-6) 

Proofs of most of the results stated in this section may be found in 
Bland [12,1968, pp. 79-90]. 
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and lead to two generalized Fourier series, one for downwash, the other 


for pressure. Let 

<f '9>w = 7 £ i f(x)g(x)dx, (3-7) 

<f '5>p = J £* /^§ f( C)g(C)d? (3-8) 

denote complex inner products with respect to the above weight functions, let 

||f|| w = (3-9) 

||f|| p = Af ,f>p (3-10) 

denote their inner product norms, and let 

Lw = {f! ll f llw<“}r (3-11) 

= {fs||f||p<-}. (3-12) 


The equivalence classes of L 2 and L 2 induced by these norms 2 are Hilbert 
spaces and the polynomials (x^l an< ^ t^n^i can be shown to be respective 
bases for them. Thus, if fsL^ and gsL 2 , then 

ir 

f = n=i f Xn Xn ' 

9 = nil g K 'Pn' 

where we adopt the abbreviated notation 

T_ 1 / 1+X 

f x n E <f 'X n >w = 7 / l^r f (*>x n (*> d *' 

= <9'*n>p = V / /j^ g(E)* n (5)d5 

for the generalized Fourier coefficients. Since these polynomials are 
complete, the only functions which are orthogonal to all of them are null; 


(3-13) 

(3-14) 

(3-15) 

(3-16 


2 The equivalence class of feL^ is the set of all functions g such that 
1 1 f — g 1 1 w = an< ^ similarly for feLp. Xn general the Fourier series of a 
function need not converge everywhere to the value of the function itself, 
particularly in the case of a flap where the downwash function is discon- 
tinuous. The notion of an equivalence class enables functions and their 
Fourier series to be treated as essentially the same. Henceforth, unless 
stated to the contrary, we shall not distinguish between different representa- 
tives of equivalence classes. 
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x.e. , 


fel/ & f Y =0Vn+f=0, 

W A n 

geLp & g^ n = 0Vn->-g=0. 

The role of airfoil polynomials in aerodynamics stems from their 
being integral transform pairs of one another: 

1 r 1 .O+x X n ( x ) ^ 

jU /i^^r 6x = 

IT -1 / 1+5 x-5 Xn(x) 

and from recognizing that the airfoil equation and the Sohngen inversion 
formula can be recast in the same form: 


(3-17) 

(3-18) 


(3-19) 

(3-20) 


P(?) 


IT _1 / 1— X X-5 


4w(x) = 1 / pm. 

7T _i v 1+5 X-5 


From the above, it follows that 

OO Q 

w = E w- v r - z 


2 w V X & PeL z P = 42 w Y & weLA. 

n=l ^n A n P n=l x n r n ^ 


dx. 

(3-21) 

dC. 

(3-22) 

■ n i|» n & weL 2 . 

(3-23) 


In other words, for steady incompressible flow, the generalized Fourier 
coefficients of the nonsingular pressure function P with respect to {^ n )i 
are equal to four times the generalized Fourier coefficients of the downwash 
function w with respect to t X n J l - i-e., 

P* n = <p ^n > p = 4 <w 'X n > w = 4 w x n * < 3 - 24 ) 

For this reason, the functions x n are called downwash polynomials and the 
functions i|; n are called pressure polynomials . 

The integral transforms appearing in (3-19) - (3-22) are central to 
the theory of airfoil polynomials. They appear to have been first discovered 
by N.I. Akhiezer in 1945 and are discussed briefly in [9,1976,Ch 2,p,133j. 
However, since these insights were apparently initially applied to the 
solution of practical aerodynamics problems by Bland in his thesis [12,1968] , 
we shall define the integral transforms 
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1 1 

Hf (x) = - / 

TT -! 

H -1 f (?) = - / 

TT _ X 


/iri 

/ 1+5 

/l+x 
/ 1-x 


fjj) 

x-5 


f (X) 

x-5 




dx 


as Bland transforms . 

The functions x n anc ^ are related to Chebyshev polynomials and can 
be shown to be constant multiples of Jacobi polynomials. 3 Their n-1 zeros 
are given by 

2i7r 

Xn ) = 0 ; x£ — —cos 2.r\~2. ' ^ = ^ * ■ • • * ri — 1 ? n 2 , 

*n<5S> = 0; 5? = cos i = n > 2, 

are antisymmetric in the sense that £? = -x?, and are interdigitated 

according to 

-i < e5_! < x? < ••• < e? < < i. 


as shown in Figure 4. 


(3-25) 


(3-26) 


( 3-.27) 
(3-28) 


(3-29) 



Figure 4. Interdigitation of zeros 

3 See, Bland [ 12 , 1968 j?p. 79-80] and Abramowitz and Stegun [13, 1964, Ch. 22]. 
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F 


The logarithmic transforms 


i LJ X n (x)log|x-c|dx = 


\ ^ n (5)log|x-5|dC = 


permit closed form integrations of a logarithmic singularity appearing 
in the kernel, and the continuous part of the kernel may be integrated 
according to the following NQ-point Jacobi-Gaussian quadrature rule: 


X? (x) + (1+log 2 ) y i (x) ' 
2 


n = 1 


X n +] (x)+x n (x) _ XnW+Xn-l (*> n > 2 

2n 2 (n-1) ' — 


(3-31) 


_ i^2 (C) + (1+log 2)*! (5) 

2 


n = 1 


»n+i<g>-»n(€> fa(€)-'Pn-i(g) 

2n 2Tn-l) ' 


(3-30) 


/ /l3 f(x)dx “ ^ (i+xNQ + i)f( X NQ +1 ), 


2NQ+1 i=l 

1 /' 1-C - . , , „ 2n 


NQ 

2 


/lif f<5>aE ■ 255a n (i-e;®* 1 * * )*!?! 6 * 1 )- 


(3-32) 

(3-33) 


If the generalized Fourier coefficients for the continuous factor 
of pressure are known. 


p = X P* n #. n (3-34) 

then the lift and pitching moment coefficients 4 * 

1 1 

C L = j- f Ap (£ ) d£ , (3-35) 

1 1 1 

C M = 2 / (£+y)Ap(£)d C (3-36) 

are given exactly by the first two Fourier coefficients alone: 

C L = /ir|'"l (?) nll p * n ^(C)dC = f P#,. (3-37) 

C M = 7 ^ / irf p * n ^n(?)<35 = J P^ 2 . (3-38) 


4 In this definition lift is positive upward, and pitching moment is positive 

counterclockwise (leading edge down) about the quarter chord and is based on 

the semichord. 
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Ex act closed form expressions similar to (3-37) and (3-38) may be 
obtained for the generalized aerodynamic force matrix. The components 
A rs of this matrix are defined as the virtual work done as the structure 
deforms in mode r against the pressure corresponding to mode s: 

1 1 

A rs E T f h r (5)Ap s (?)dC; r, s = 1,2,... (3-39) 

Z —i 

where {h r }r=i are structural basis functions, commonly selected in practice 
as the in vacuo vibrational eigenfunctions. Expand the structural basis 
functions and their corresponding pressure functions in the pressure 
polynomials 


hr(£>) — < ^ 1 r , ^n > p ^n^^ — n=i ^-nr ^n^£) ' r “ 1,2,... 

A Ps(5) = X <Ps^n > p *n«) = P »s 'MO; s = 1,2,... 


and obtain 


1 00 


A rs 2 n=l < ^ 1 r ,l ^n > p <p s ,l ^n > p 2 n=l ^nr Pns? r ' S ^ ^ 1 


(3-40) 

(3-41) 


(3-42) 


The generalized aerodynamic force matrix may be considered to be the single 
most important aerodynamic quantity in aeroelasticity because it completely 
determines the aeroelastic coupling. 
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§4. Analytical and numerical properties of Bland 's equation 

In this section we discuss some of the theoretical aspects of the 
collocation method developed by Bland for the solution of (2-57) . In 
particular it will be shown that the apparent numerical convergence 
observed by Bland and ourselves can be established theoretically. In 
order to do this it is necessary to take a somewhat different viewpoint 
toward collocation than has been done in the past. We regard, as did 
Bland, collocation as only a secondary numerical method arising out of 
somewhat different procedures. Bland, in his thesis, proposed a general- 
ized least squares method for the solution of (2-57) , which was shown 
under certain conditions (to be discussed later) to be equivalent to 
collocation. In the same spirit we approach Bland's equation by a 
projection method analogous to Galerkin's method for equations of the 
second kind. Here again it can be shown that this technique is numerically 
equivalent to collocation. Reversing the argument, we prefer to regard 
the Galerkin method as primary, since it will be shown that the theory 
appears to emerge best from this viewpoint. 
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4.1 Collocation 

For completeness and ease of reference Bland's collocation method 
is now reviewed. Since the space Lp is spanned by the pressure polynomials 

oo 

{ ^ n ) 1 F p can k e expanded as 

P = J, <P^n> P V 

because of this one seeks an approximate solution P N of the form 

N 

P N = n ii a n fer (4-1) 

where a n ; n = are to be determined. Letting 

1 /T^i 

r N (x) = w (x) - / K(x-5)P N (5)d5, (4-2) 

and setting r N (x) equal to zero at the zeros of x^j+i ( x ) gives N linear 
equations in the N unknowns (a n ]q. Solving for these determines the approxi- 
mation P N to the pressure coefficient P. This is Bland's collocation method. 
The numerical examples given by Bland and by us appear to indicate that P N 
converges to P. 

To establish this it is convenient to regard (2-57) as an operator 
equation acting between the two spaces Lp and L £ defined in Section 3. 

Making use of the properties of the airfoil polynomials given by (3-15) , 

(3-16) , (3-30) and (3-31) enables us to show that this equation has interesting 

analytic properties which permit a detailed analysis of the above numerical 
method. 
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4.2, Analytical properties 

To set the stage for our theoretical treatment of (2-57) we make 
some simple preliminary observations. From (2-50) the kernel can be 
written as 

K(x-£) = log | x-5 | +K c (x-£), (4-3) 


where K^>(x-C) is the bounded part. Multiplying both sides of (2-57) 
4 

by — allows us to rewrite it as 

p 


•l-C 


w(x) = l, ✓ 


K(x-£)P(£)d£, 


(4-4) 


4 4 

where w(x) = (— )w(x) and K(x-£) = (— )K(x-£). This transformation obviously 

p p 

has no effect on the solutions of (2-57) . The key observation, and one 
that was made by Bland in his thesis is that (4-4) is actually equivalent 
to an equation of the second kind. To see this decompose K(x-£) as 


K(x-£) 


1 

tt(x-£) 


+ K sc (x-C) . 


(4-5) 


This first term in (4-5) leads to the integral operator 


H P (x) 


1 

TT 


/ A -e p(£) 

-1 / 1+E x-£ 


d£ 


which is just the Bland transform (3-25) . From (3-19) and (3-20) it is 
seen that H maps the orthonormal basis {^ n }l onto the orthonormal basis 
(Xn^l' and thus can be extended to a bounded invertible operator from 
to L^. If one then operates on both sides of (4-4) by H -1 it takes 
the form of an equation of the second kind in L^. This leads one to suspect 

xr 

that the well developed numerical methods for such equations would be useful 
in solving (2-57) which, as we shall see, is the case. 

To proceed further, it is necessary to develop some of the analytical 
properties of (4-4) . For this we digress slightly and state some definitions 
and theorems from functional analysis. 


Notation : Let Hj and H 2 be Hilbert spaces. The inner products on Hi,i = 1,2 

will be denoted by < , >^ and the corresponding norm by || || j_. The set of 

bounded linear operators from to H 2 will be denoted by [H^,^]. 


Definition 4.1. Let Hj and H 2 be Hilbert spaces. Let He[Hi,H 2 ]. If 
||Hx|| 2 = || x Hi for all xeH, we say that H is an isometry . In addition if 
H has a bounded inverse we say that H is unitary, 


-31- 


Theorem 4.1. Let H , be a Hilbert space with a complete orthonormal 
basis { ^ n } i . A necessary and suff ic ient condition that He[Hi,H?] be unitary 

CO 

is that the s et be a complete orthonormal basis in H? . 

Proof. See [81] . 

Definition 4.2. Let Hi and H 2 be Hilbert spaces. Let Kg[Hi,H 2 ]. We say 
that K is compact if for every bounded sequence {x n }i in Hi the sequence 
{ Kx n ) i has a convergent subsequence. 

Theorem 4.2. (Properties of compact operators.) Let Kj£ [Hi ,H 2 3 ? i = 1,2 
be compact. Then 

(i) Ki + K? is compact. 

(ii) If Te[H],H?] then T K-j ; i = 1,2 are compact. 

(iii) If C is a complex number then CKj ; i = l y 2 are compact. 

(iv) If Ki has finite dimensional range then Ki is compact. 

(v) Let |1 || denote the operator norm on [Hi ,H?] . Let {K n } 7 be a sequence 

of compact operators from Hi to H 2 such that lim ||Kn-K[f = 0 where Ke[Hi,H 2 ]. 

n-*” 

Then K is compact. 

oo 

(iv) It follows from (iv) that if {K n )i is a sequence of operators each 
having finite dimensional range, and lim ||K-K n || = 0 then K is compact. 

Proof. See [81] . 

Definition 4.3. Let He[Hi,H 2 ]. The adjoint of H denoted by H* is the 
bounded linear operator defined by 

<H*x,y>! = < x ,Hy> 2 , (4-6) 

where x and y are arbitrary vectors in H 2 and Hi respectively. 

Theorem 4.3. Let He[Hi,H 2 ]. Then H has a unique adjoint H*e[H 2 ,H~|]. 

Proof. See [81] . 

Theorem 4.4. Let H be a unitary operator from Hi to H? . Then H* = H -1 . 

(H - * 1 denotes the inverse of H.) 

Proof. Note first of all that since H is an isometry it preserves inner 
products; i.e. , <Hx,Hy> 2 = <x,y>i for all (x,y)eHi- In fact ||H(x+y)|| 2 

= II (x+y) II i 


But 
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||H(x+y) || 2 = <Hx,Hx >2 + 2 Re <Hx,Hy> 2 + <Hy,Hy> 2 

= <x,x>! + 2 Re <Hx,Hy> 2 + <y,y>i • (4-7) 


Similarly 

|| (x+y) || ^ = <x,x>j + 2 Re <x,y>j + <y,y>i- (4-8) 

Equating (4-7) and (4-8) along with replacing x+y by x+iy gives the result. 

NOW 

<Hx,Hy> 2 = < H*Hx f y>! = <x,y>!, (4-9) 

by the definition of adjoint and the above observation. Since (4-9) holds 
for all yeHi we get that 

H*Hx = x, VxeHj. (4-10) 

Let yeH]^ Then, since H*" 1 exists, y = Hx for some xeHi . Therefore 
y = Hx = H (H*Hx) = (HH*)Hx = HH*y. Thus y = HH*y for all yeH 2 , and this 
with (4-10) shows that H* is both a left and right inverse of H so that 
H -1 = H*. 

00 

Theorem 4.5. Let Hi and H? be Hilbert spaces. Let Te[Hi,H?]. Let {^ n )i 
be a complete orthonormal basis for . Assume that 

INn||| < (4-11) 

Then T is compact. 

(Note: Operators satisfying (4-11) are said to have finite double or 

Hilbert- Schmidt norm. See [82] for a discussion of their properties.) 

Proof. From Theorem 4.2 it suffices to show that there exists a sequence 

of operators {T n }j, each having finite range, such that lim 1 1 T-T n 1 1 = 0. 

nr+°° 

Let U n = span and define Q n by 

n 

Q n x = k Z l <x,^ k > 1 i/; k , xeHi. (4-12) 

Q n is the operator of orthogonal projection onto U n . Let T n = TQ n . From 
(4-12) it is seen that 

n 

T n x = k i L <x /4>k > l T ^k' xeH l- 


(4-13) 
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Letting U n = span {Ti^}] 1 , we see that T n has finite range. Now 

1 1 Tx-T n x 1 1 2 = l! k= J +1 <x »*k>l 1 1 2 1 k= | +1 l<X,l(» k >x| l| Tt|/ k || 2 • (4-14) 

By the Cauchy- Schwarz inequality (4-14) becomes 


I Tx-T n x | 


2 1 


( k=n+l 


| < x,i|) k >i 


l 2 ) 2 ( 


, l 
k=n+i 


INfellS) 


2 — 


1 ( k=n+l 


IN k ll|> 


(4-15) 


From (4-15) it is seen that 

1. 

||T-T n || <. ( k J +1 ||T^ k || |) 2 . (4-16) 

Since Z ll^kHI converges the right hand side of (4-16) can be made 
arbitrarily small for n large. Thus lim 1 1 T— T n 1 1 = 0 and the theorem is 
proved . 

The results of the above theorems are now applied to (2-57) (or 
equivalently (4-4)). As was stated in §4.1 we wish to regard (2-57) as 
an operator equation acting between L p and L^. Using the equivalent form 
(4-4) the kernel decomposition gives rise to three operators H, and K 2 
defined by 


H P(x) 



1-C P(5) 
1+E 


d€. 


(4-17) 


K x P (x) 


ik 

7T? 2 " 


i-e 


-1 S 1+5 


log|x-C|P(5)dC, 


(4-18) 


and 


1 l-£ 

K 2 P(x) = / y K c (x-C)P(^)d^ / 


(4-19) 


where the integral in (4-17) is taken in the sense of a Cauchy principal 
value. The basic properties of these operators are summarized in Theorem 4.6 
below. 

Theorem 4.6. H, Ki and K? define bounded linear operators from Lp to l£. 

In addition H is unitary and Ki and K? are compact - 

Proof. As stated above H is just the Bland transform (3-25) . Since 

OO OO 

H0 k = Xk where and (Xk^l are defined by (3-1) and (3-2) respectively, 

H can be uniquely extended as a bounded operator from Lp to L^ by 

HP(x)= k | 1 <P,^ k > p H^ k (x) = k | x <Pf^k> p X k - 


Since (Xfc} is a complete orthonormal basis in L 2 it follows from Theorem 4.1 
that H is unitary. 

From (3-31) we see that 

X?(x) + (1+log 2) xi (x) _ 1 

2 / n — ± 

Xn+1 + Xn _ (Xn+Xn-j) n S- 2 

, 2n 2 (n-1) “ J 

Since is well defined on the basis elements {^ n )i it can be extended 
as an operator to all of by the formula 


Ki^n = 


ik 

tTb 2 " 


KiP(x) = k | 1 <P/^ k > p (4-20) 

provided that the sum in (4-20) converges in the L 2 ^ norm. To verify this 
observe that 


l|Kl*nllw = <K l*n' K l*n > w 


2+2 log 2 + (log 2) 2 
k 2 I 4 

7t 3 4 (n 2 -n+l) 

2n z (n-1) 2 ' 


n = 1 
n > 2 


Thus (4-20) converges, and since | j 1 1 ^ is of order the series 
n ?^ || K i^nllw converges and thus by Theorem 4.5 ^ is compact. 

Now 


< 


|k 2 P(x) I 1 / / |Kc(x-C) I |p(5)|d£ 

/£§ ME)l 2 dElV /if iKcte-Ol^El 2 


by the Cauchy-Schwarz inequality. Thus 

I|k 2 p|Iw 1 W p \\p l\ /ill [/ /'ill |Kc(x-C)| 2 dC]dx, (4-21) 

where the inte gral in (4-21) exists since Kc(x-£) is bounded and the 

/ 1+x 

— — is integrable on [-1,1]. From this it is seen that 
||k 2 p|| w 1_ C||p|| and so K 2 is bounded. To establish compactness we again 
resort to Theorem 4.5 by showing that ll K 2^ n Hw conver 9 es - For 

observe that 

Ji l^kl 2 = k Ii |/ /ii|Kc^-?)^k(5)d?| 2 . 


(4-22) 
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But 


where n x (£) 



K c (x-C)^ k (C)d^ = it <n x ^k > p^ 


(4-23) 


KcCx-C). Thus the right hand sum in (4-22) can be written as 


" 2 jJj l <n x'^k>pl 2 - (4-24) 

By Parseval's Theorem thi s is equal to tt / J y— |r) x (£)| 2 d£. Multiplying 

/ i+x """ 1 ^ 

both sides of (4-22) by J and integrating with respect to x gives 



11 * 2 * 11 * 



[//|^||Kc(x-?) | 2 d?] . 


(4-25) 


/ 1+x 

The integral in (4-25) is finite since K^(x-C) is bounded and J is 
integrable. Thus K 2 is compact since it has finite double norm. 

From the results of Theorems 4.5 and 4.6, (4-4) can be written in 
operator form as 


(H + K) P = w (4-26) 

where H + K e [Lp,L^] and K = + K 2 is compact. A solution to (4-4) will 

now mean an element PsLp solving (4-26) . In general such a solution will 
satisfy (2-57) or (4-4) almost everywhere. If the solution is sufficiently 
smooth it will satisfy (2-57) in the usual pointwise sense. Since physically 
one is interested in weighted integrals of P such generalized solutions 
are perfectly reasonable. 

From Theorem 4.4 we know that H has an inverse. Applying to both 

sides of (4-26) gives the equivalent equation 

P + H" 1 K P = H" 1 w. (4-27) 


From Theorem 4.5 H” 1 K is a compact operator and (4-27) has the form 

(I + L) P = H — 1 w (4-28) 

where I is the identity operator on Lp and L = H _1 K. Equation (4-28) is 
now in the standard form of equations of the second kind [83] . From 
this it is possible to obtain solvability theorems by appealing to the 
Fredholm theory. We state only one such theorem. 

Theorem 4.7. Bland 1 s integral equation has a solution in Lp iff -1 is not 
an eigenvalue of H" 1 K. 
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Proof. This is just the usual Fredholm alternative for operators of the 
form I + L where L is compact. See [84] for more details. 

Throughout the rest of our discussion it will be assumed that the 
solvability condition in Theorem 4.7 is satisfied, and we now proceed to 
a discussion of the numerical solution of (4-26) . 
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4.3. G a lerki^s method for a class of operator equations 

Motivated by the results in 4.2 we now consider the numerical solution 
of the class of equations 

TP = w (4-29) 

where Te[Hi,H 2 ]r H^; i = 1,2 are Hilbert spaces, T = H + K where H is unitary 
and K is compact. From Theorem 4.4, (4-29) can be written in the equivalent 
form 


(I + H*K)P = H*w. 


(4-30) 


(4-30) is an equation of the second kind and since there are many well 
understood methods for solving such equations [84] we consider the 
possibility of using them to find numerical algorithms for (4-30) . For 
our purposes Galerkin's method appears to make best use of the theoretical 
structure of (4-30) and we now give a brief discussion of this popular 
technique. 

oo 

Assume that {^ n }i is a complete orthonormal basis for Hj and look 
for approximate solutions of (4-30) of the form 

P N = n Ij a n i|» n , (4-31) 

where the a n 's are constants to be determined. Since in general P N will 
not solve (4-30) exactly we consider the residual Rn given by 

% = Pn + H*K P n - H*w. (4-32) 

If Pn were the true solution of (4-30) then R n = 0. However, this will 
not be the case in general and we try to pick Pn so as to make small. 
Galerkin's method attempts to do this by making R N orthogonal to ^ n ; 
n = 1,2,...,N. That is we require 


^ Rn *^ n > l — 0 1 n 1,2,... ,N. 


(4-33) 


Substitution of (4-32) into (4-33) gives the following set of linear equations 
to determine the a n 1 s : 


N 

a n + a m = < H*w,iJ>n > i/ n = 1*2, ...,N, 

where the orthogonality of the ^ n 's ^ as been used in the derivation of 


(4-34) 


N 


(4-34). If the equations in (4-34) have a unique solution then a n ; n = 1,2,..., 
can be determined and the approximation Pn is well defined. The following 
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theorems, taken from Atkinson [84] justify this procedure. First we 
recast Galerkin' s method into a slightly more abstract form. 

As in Theorem 4.5 let U N = span {i/ ; n ]q and let Qn be the operator of 
orthogonal projection onto U N - The Galerkin equations (4-33) and (4-34) 
are equivalent to the operator equation [1976,1] 

Qn r n = °/ p N£%- (4-35) 

Or, written out in full 

Qn( I+H*K)P n = Q n H*w, P N eU N . (4-36) 

Assume that (4-30) has a unique solution. Then from the Fredholm 
alternative [84] the operator I + H*K has a bounded inverse since 
H*K is compact. 

Theorem 4.8. Let % = H*K. Then ^im ||k-K n || = 0. 

Proof. See Atkinson [84] . 

Theorem. 4.9. Let Kn be as in Theorem 4.8 and assume that N is large enough 
so that ||k-K n || < I I i+h*k | | " T ^ en ^he operator (I+K N ) _1 exists, is bounded and 


II (I+Kn> Ml 


(I+H*K) 


-1 


- 1- I+H*K 


K k N 1 1 


(4-37) 


From this it follows that the Galerkin equations (4-36) have a unique solution 
and 

l|p-%lll 1 II (I+Kn) _1 |I I|p-QnP|Ii- (4-38) 

(4-38) implies that P N converges in norm to P the solution of (4-30) . 


Proof. See Atkinson pp. 51-52 for details of (4-37) and (4-38) . The con- 

oo 

vergence follows from (4-38) and^the fact that { 1^)1 is complete in Hi so 
that || P— QnP II 1 = ^ N Ji H p '^ n >:L I ^ which converges to zero. 

Note that (4-37) and (4-38) constitute an error estimate for the approxi- 
mate solution P N . 

Theorem 4.9 is the basic convergence result that we will use in discussing 
Bland's method. Note that it has the immediate consequence of showing that 
the application of Galerkin 1 s method to Bland's equation (2-57) is a con- 
vergent numerical method. That this is true follows from the fact that 
the equivalent operator version (4-26) satisfies all of the hypotheses to 
prove Theorem 4.9. 
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As pointed out in Atkinson there exists a "dual" to Theorem 4.9. 

That is, if one can establish the existence of a unique solution to the 
Galerkin equations for some sufficiently large N, then it follows that one 
can prove the existence of a solution to (4-30). From a practical point 
of view this is important, since we have the numerical information available 
from the TWODI program, and thus the solvability of the numerical problem 
can be used to infer the existence of solutions to the original problem. 

Although we have established that Galerkin* s method is theoretically 
a reasonable procedure to use numerically, it does initially appear to have 
several drawbacks, the most important of which is the necessity of performing 
the complicated integrations needed to evaluate the inner products in (4-34) . 
Since, practically these integrals must be done numerically it is important 
to examine the effect of this on the Galerkin equations. Surprisingly, due 
to the structure of T, particularly the unitar ity of H, considerable simpli- 
fication results, and as will be shown, under appropriate conditions Galerkin* s 
method becomes equivalent to collocation. 
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4.4. The relation between Galerkin's method and collocation 

We now proceed to examine the above mentioned relation between collo- 
cation and Galerkin's method. Recall from (4-34) that the Galerkin approxi- 
mation for P is obtained by solving the linear equations 
N 

a n + m S 1 <H * K ^m'^n > l ^ = <H*w,* n > lf n = 1/2,..., N. 

Using the definition of adjoint, these equations become 
N 

a n + hi— ^ < ^m'^n > 2 ^m = ^ — 1/2, (4—39) 

oo oo 

Since H is unitary { Hi|; n } ^ = ^ X n } i is a complete orthonormal basis 

for H 2 - Using this and the fact that < X n rX m > 2 = ^mn (4-39) becomes 

N N 

mSj f <H ^ m , X n > 2 + < K0 m ,x n > 2 } a m = <w, Xn >2 = m Ij < iH+K) ii» m ,X n > 2 a m 

= < W /Xn > 2 r ^ = 1/2, 

Or, since H+K = T , 

N 

< Ti|; m ,x n > 2 a m = < w 'X n > 2 ? n = 1/2, • . . ,N. (4—40) 

Note that (4-40) is formulated directly in terms of the original equation 
of the first kind TP = w and can be regarded as a projection method in its 
own right. To see this, again look for an approximate solution to (4-29) 
in the form P^ = n f 1 a n i|; n . Using the same argument as for Galerkin's method 
leads us to consider the residual 

r N = TPjj-w. (4-41) 

In order to make r^ "small" we pick a n ; n = 1,2,...,N to satisfy the or- 
thogonality condition 1 

< r^ ,Xn > 2 = n = 1,2,...,N. (4—42) 

Writing (4-42) out in full gives (4-40). If one knows (Xn^? explicitly 
then numerically it is more efficient to use (4-42) than (4-40) . 

Since our main interest is in Bland's equation (2-57) we now specialize 
(4-40) to the case where H^; i = 1,2 are taken to be Hilbert spaces of 

^his is sometimes called the method of weighted residuals [85] . 
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functions on the interval [a,b] and the operator T is an integral operator 
of the form 

TP(x) = /g K(x,e)P(£)d€, (4-43) 

where the integral in (4-43) is generally taken in the sense of a Cauchy 
principal value. Let 

L^x) = Tifr m (x). (4-44) 

Then 

<T *m»Xn > 2 = <L m'Xn > 2- (4 " 45) 

Since the inner product in (4-45) is an integral, we assume that it can 
be approximated by a quadrature rule having N nodes, denoted by x k ; 
k = 1,2,...,N and corresponding weights W k ; k = 1,2,. . . ,N. Thus (4-45) 
can be approximated by the finite sum 

N 

^ x k^n ( x k) • 

A similar approximation applied to < w,x n > 2 gives 

N 

<w 'Xn > 2 ~ k i x W k w(x k )x n (x k ) . 

Thus the numerical solution to the Galerkin equations (4-40) is obtained 


by solving 

WkLm^kJXn^ )) = k Ii W k w(x k } X n < x k> - n = (4-46) 

Define vectors and matrices by 

a = { a n ) ; w = { w(x k ) } (4-47) 

£ = tX n (x k ) 1 , L = [Lm^k) 1 # W = diag [W k ] ; m,n,k = 1,2, ... ,N. (4-48) 

Using (4-47) and (4-48) , (4-46) takes the form 

(^WjLa^ = (^W)w. (4-49) 

If it is assumed that the matrix x_ is nonsingular then is nonsingular 
and (4-49) is equivalent to 

La = w. (4-50) 


Referring back to Section 4.1 it is easily seen that for Bland's equation 
(4-50) are just the collocation equations resulting from (4-2) . Thus it 
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is seen that if the quadrature errors are ignored in the solution of 
the Galerkin equations then they will give exactly the same numerical 
solution as collocation, provided that the matrix x is invertible, as 
will be established below. Thus from a numerical standpoint it makes 
little difference whether one uses collocation or Galerkin* s method in 
the solution of (2-57) — the result will be the same numerical approxi- 
mation to P. However, from a theoretical point of view Galerkin* s method 
appears to be more desirable, since as we have seen, it gives convergence 
along with computable error bounds. 

To complete our equivalence proof of collocation and Galerkin* s method 
we now establish the nonsingularity of x* 

Definition 4.4. Let {f n (x)}; n = 1,2,. . . ,N be N functions defined on the 
interval [a,b] . We say that {fn<x>}? are unisolvent if the matrix [fn( x k)l? 
n,k = 1 , 2 , . . . , N is nonsingular for every set of distinct points {x^.}^ in [a,b] . 

Theorem 4.10. Let (p n (x) be a basis for the polynomials of degree <_ N-l 
on [a,b] . Then (p n (x)}^ are unisolvent. 

N-l , 

Proof. Let tt = [p n (xk)l; n,k = 1,2,...,N. Since p(x) = , p n k xK it follows 

X— 0 

that tt = PV where P = [p n k^ and v is the Vandemonde matrix given by 



i 

1 . . . 

. 1 




x i 

x 2 • . . 

• X N 



V(Xi ,x 2 , - . .x N ) 



o 


• 


X? 

x£ . . . 

* X N 




xN- 

1 X?- 1 - . . 

x N ‘ 
• X N 

-1 


Thus det tt = det P det V. But det 

P is nonzero 

since 

{ p n (x) is a basis 

det V = (-1) N II 

(x L - Xj) ? 

o. 



i< j 




if j=lf 

• . . , 

N 




since the points x^ are distinct. 

Thus det tt ? 

0 and 

tt is nonsingular. 


Corollary 1. Let {p n (x))i J be polynomials of degree <_ N-l orthogonal with 
respect to some inner product on the set of functions on [a,b] . Then they 
are unisolvent. 
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Proof. The orthogonality implies that {pnCx)}? are linearly independent 
and thus a basis for the polynomials of degree <_ N-l on [a,b] . The corollary 
now follows from the theorem. 

Corollary 2. Let (x n }? ^e the first N downwash polynomials as defined in 
Section 3 . Then { x n )? are unisolvent. 

Proof. Since (x n ^ are orthogonal on [-1,1] with respect to the inner 
product < , > w by corollary 2 they are unisolvent. 

Corollary 3. Let (x n J? be as above and let {x^}^ be the zeros of X n +i ' 
then the matrix ^ = {Xn( x k)} is nonsingular. 

Proof. From (3-24) we see that X^+i ^ as N distinct zeros on [-1,1]. Since 
{ Xn }? are unisolvent the result follows. 

From Theorem 4.10 and the above discussion we arrive at our main equiva- 
lence result for collocation and Galerkin's method for Bland's equation. 

Theorem 4.11. Let P N be the approximate solution to (2-57) given by using 
Galer kin's method^ based on the pressure po ly nomials as a complete ortho- 
normal basis for L|^. Then if the inner products in (4-40) are evaluated 
using the Ja cob i-Gaussia n quadrature rule (3-29) , the resulting numerical 
approximation^ to^ (2- 57) is the same as one obtains using collocation with 
the same basis and collocating at the zeros of x N+1 (x) . 

Proof. From (4-49) it suffices to prove invertibility of x_ which was done 
in Theorem 4.10. 

The importance of Theorem 4.11 is that it enables us to regard collocation 
as numerically equivalent to Galerkin's method. As we have seen Galerkin's 
method is convergent, and since neglecting the quadrature errors in the 
evaluation of the inner product in (4-40) gives Bland's collocation equations, 
we can conclude that this method is convergent. We summarize this, our main 
result, as Theorem 4.12. 

Theorem 4,12. Let P^ be the numerical approximation to the solution of 
(2-57) as described in Section 4.1. Let P^ be the Galer kin approximation 
to P as described in Section 4.3. Let Pg^ be t he approximation to pG obtained 
by evalua t ing the inne r products b y the Jac obi-Gaussian rule (3-29) . Then 
P N = P N and thus neglecting these quadrature errors pG converges in the norm 
of Lp to P and thus so does Pg The error in this approximation is given by 
(4-38) and so 
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/' 


1 

l|PN- p Hp - C N H p -2N P Hp = C N< n J +1 l <p ^n>p| 2 ) 2 - (4-51) 

Note that (4-51) shows that the rate of convergence of P§ to P depends 
on the smoothness of P and thus requires a knowledge of the behavior of 
the generalized Fourier coefficients of functions expanded in airfoil poly- 
nomials. In addition it requires estimates of the smoothness of the solu- 
tions P. We expect to pursue these points in future work. 

As we stated at the beginning of this section Bland's starting point 
for the solution of (2-57) was a least squares procedure which he claimed 
was equivalent to collocation. His proof of this fact [12] resided in 
the assumption that the collocation matrix (L n (xk)} was nonsingular. As 
a further consequence of the Galerkin theory we establish the validity of 
this proposition. 

Let G denote the matrix [< ^ X n > 2 ^ given in (4-40) . From the proper- 
ties of Jacobi-Gaussian quadrature < Ti/j m ,Xn > 2 = 9nm + e nm where [g^] is the 
matrix x_WL. as defined in (4-49) and [e^] = E is the matrix of quadrature 
errors. Thus 

G = x_WL + E. (4-52) 

By Theorem 4.10 

L= (XW) _1 (G-E). (4-53) 

From (4-53) it is seen that L has an inverse if and only if G-E does. From 
Theorem 4.9 we know that for all sufficiently large N, G has an inverse. 

Thus G-E = G(I-G _1 E) and it follows from Banach's lemma [85] that 
I-G _1 E has an inverse provided that ||g“ 1 e|| < 1 where || || is any matrix 

norm on the set of NxN complex matrices. Since ||g _1 e|| <_ | [ G — 1 1 1 ||e|| it 

suffices to show that ||e|| < n - 1 r r r . Under the assumption that the quad- 

11^ II 

rature errors can be made arbitrarily small if N is large enough we can 
pick N so that G” 1 exists and ||e|| < j| . Thus we conclude that for N 

large enough the collocation matrix is nonsingular and consequently Bland's 
observation that collocation is equivalent to least squares is valid. 

Theorem 4.13 states this result. 

Theorem 4.13. Provided that N is sufficiently large Bland's l east squares 
method, his collocation method and our Gal erkin method all give the same 
numerical approximation to the solution of (2-57) . 
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4.5. Conditioning of the collocation matrix 

As a final application of the Galerkin approach to (2-57) we offer 
a brief discussion of the conditioning of the collocation matrix L. In 
[12,1968] and [5,1970], Bland remarks on the fact that L is well conditioned 
without offering any proof. The accuracy of our own numerical results also 
supports this observation. It also appears to be part of the folklore of 
singular integral equations of the first kind that the strong Cauchy singu- 
larity leads to numerical methods which are well conditioned [86,1971] and 
[87,1975]. Although this is intuitively reasonable we are unaware of any 
mathematical proof of this fact. For equations of the second kind there 
exist fairly complete results on conditioning. A summary of these may be 
found in Atkinson [84,1976]. Since we know that (2-57) is equivalent to 
an equation of the second kind these results should be of use here. That 
this is the case is demonstrated below. 

Definition 4.5. Let A be an NxN complex matrix. Let || || be a matrix norm 

on C N . The condition number of A relative to || || is given by 

C (A) = II A II II A“ 1 || . (4-54) 

A matrix is said to be well conditioned if C(A) is of order 1 and 
poorly conditioned if C (A) is large [84,1976]. it follows immediately from 
(4-54) that if A and B are matrices then 

C (AB) < C (A) C (B) (4-55) 


and 


C (A -1 ) = C(A) . 

Since L = (x.W) _1 (G-E) = (^W) -1 G -1 (I-G _1 E) , (4-55) and (4-56) give 

C (L) <_ C(^)C(W)C(G)C(I-G -1 E) . 


(4-56) 


(4-57) 


Since W is a diagonal matrix and one can show that C(W) <_ c 


max 


Wk 


Wk 


where 


c is a constant independent of W. Since W^; k = 1,2,...,N are the quad- 
rature weights C(W) - cN. Now for N large I-G _1 E - I so that C(I-G~ l E) 

- C ( I ) - 0(1). Thus C (L) 'C , NC(x)C(G) . From Atkinson’s results [1976,1] 
one can show that ||g|| <_c M so that C (L) ~ c' n NC (x ) , which gives a reasonable 
estimate of the conditioning of L. It appears then, if N is not too large, 
that C (L) is well conditioned. A more complete analysis will have to await 
future work. 
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4.6. Convergence of integrated aerodynamic forces 

Although the above theory gives convergence of P N to P in the norm 
of Lp it is important to note that this generalized convergence of P N gives 
strict convergence of integrated aerodynamic forces to their true values. 

To see this we define 

Ap N (0 = /j~ P N (C) (4-58) 

as our approximation to Ap(£). Let f(£) be a real valued weighting function 
as in equations (3-35), (3-36), (3-39), etc., and Let F represent an inte- 

grated force, 

F = / Ap(e)f(5)d£. (4-59) 

-1 

We take as our approximation to F 

1 

F n = / Ap N U)f(£)d£. (4-60) 

By definition of the inner product on Ip, 

F n = tt <P N ,f> p . (4-61) 


Theorem 4.14. ^im F N = F. 

Proof. 

|f-F N | = * |<P,f>p - <P N ,f>p| = 71 |<P-P N ,f>p|. (4-62) 

By the Cauchy-Schwarz inequality 

|<P-P N ,f.> p | <_ || P-PN Up I! f II p- (4-63) 

Thus the result follows from Theorem 4.9. 
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§5, Organization of the computer program 


The computer program developed during this work is a user-oriented , 
working, pilot version written in extended FORTRAN IV. To accomplish 
this, several helpful guidelines were followed during its planning and 
development : 

(1) Input and output quantities were established first. These 
were selected primarily according to the general needs of the 
intended user community; 

(2) A hierarchial, modularized system of executive and compu- 
tational subroutines together with a complete list of all 
working equations was established and iterated once prior to 
actual coding; 

(3) Coding was performed with visibility as the primary cri- 
terion, facilitating changes; 

(4) All computational subroutines were carefully checked against 
independent calculations, using exact closed form special cases 
whenever possible to verify accuracy as well as correctness; 

(5) Due to the go-no go nature of unsteady kernel function pro- 
grams, computational correctness was deemed more important 
than algorithmic efficiency; 

(6) The organization of the computer program was structured to 
facilitate correctness of the pilot version and future modifi- 
cations. 

The calling hierarchy of the computer program is shown in Figure 5, 
with arrows indicating the direction of call. The main program, named 
TWODI for two dimensional, is strictly an executive program. TWODI directs 
two supervisory subroutines named PREP and SOLVE, and one interim subroutine 
named CHECK. If called by TWODI, CHECK systematically calls each of the 
major computational subroutines and checks their current computational 
results against the most accurately known results, which are stored within 
CHECK. CHECK is primarily a developmental tool, secondarily a diagnostic 
tool, and does not appear in released versions of TWODI . 

PREP is a supervisory subroutine, whose purpose is to prepare input 
data in a form acceptable to the solution process which is subsequently 
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implemented by SOLVE. PREP reads and prints back all input data. A title, 
supplied as input data by the user, is centered by subroutine CENTER, and 
the time and date, centered title and page number are printed at the top 
of each page using subroutine PAGE. PREP tests all input data except the 
flow parameters (M,k,nn and c w ) . Data lying outside their acceptable ranges 
result in the problem being deleted by PREP, with an explanatory comment. 

The input data specifies mode shapes at discrete points that are selected 
by the user. PREP collocates a linear combination of airfoil polynomials 
through these discrete points and stores the resulting Fourier coefficients 
for later use by SOLVE. Values of airfoil polynomials are calculated by 
subroutine AFP using recursion formulas, and collocation is performed using 
subroutine CSIMAL to solve complex (or real) simultaneous algebraic equations. 
Once the input data have been read, checked, printed back and prepared for 
use by SOLVE, PREP returns control of the program to TWODI, and TWODI calls 
SOLVE. 

SOLVE is a supervisory subroutine whose purpose is to solve the integral 
equation for each M,k,nn and c w case > and to compute and print the airloads 
for as many downwash modes as were given by input data. Each problem may 
have numerous cases of M,k,nH and c w . These are tested individually by 
SOLVE, and erroneous cases are deleted without affecting the others. The 
eigenvalues appearing in equation (2-54) are calculated by subroutine PICARD 
using Picard iteration. To avoid repetition, these eigenvalues are calcu- 
lated at the outset and stored for later use. Subroutine WASH calculates 
the matrix of downwashes at the appropriate collocation points. Subroutine 
BMN calculates the collocation matrix using closed form integrations for 
the Cauchy and logarithmically singular parts of the kernel, and uses Jacobi- 
Gaussian quadrature for the bounded part of the kernel. The quadrature 
points and collocation points are interdigitated according to equation (3-26) 
and the continuous part of the kernel is calculated by subroutine KC. The 
infinite series for F and F* appearing in the kernel are calculated by sub- 
routine SUM. If convergence of these series requires eigenvalues beyond 
those already stored, SUM calls PICARD as necessary. After the downwash 
and collocation matrices have been calculated, SOLVE calls CSIMAL and the 
generalized Fourier coefficients of the pressure for all downwash modes 
are thereby determined. SOLVE then calls LOADS which calculates and prints 
the particular combination of pressure, section coefficients and generalized 
forces as stipulated by the input data. LOADS calls AFP if pressures are 
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required. Output is supplied in redundant real-imaginary-magnitude-phase 
angle format using subroutine ATANC. After SOLVE has completed one case 
of M,k,riH and c w f the remaining cases are solved in succession. Upon com- 
pletion control returns to TWODI. TWODI then calls PREP for another problem, 
the solution to which is computed by SOLVE, until eventually all problems 
are solved. 

In its present form, the TWODI program works, it is believed to be 
free of error, and is convenient to use. Complete instructions for its 
use are supplied in Section 7. 
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§6. Some numerical considerations 


This Section discusses some of the numerical considerations made in 
computing Bland's kernel as given by equations (2-50) to (2-55). All compu- 
tations involve elementary operations using standard FORTRAN functions 
except the determination of the eigenvalues of equation (2-54) and the 
summation of the infinite series of F given by equation (2-51) , and its 
derivative F'. Since these computations, especially the latter, are the 
most difficult part of computing the kernel, they are discussed in detail. 

The positive solutions { X n } n=1 of the transcendental equation 

tan A + yX = 0 (6-1) 

are depicted in Figure 6(a) as the projections onto the X-axis of the inter- 
sections of a straight line with branches of the tangent function. Equation 
(6-1) occurs elsewhere in mechanics 1 with both positive and negative values 
of the parameter y being physically meaningful. In the present study, only 
nonnegative values of y are meaningful but the solution algorithm we use 
is equally valid for all real values of yz 

-oo < Y < oo. (6-2) 

The eigenvalues are well separated and satisfy the inequalities 

ir(n-j) < A n < 7T(n+y); n = 1,2,... (6-3) 

Figure 6(b) depicts a natural iteration scheme. Because of the behavior 
of the derivative of the inverse tangent function, if A^ k ) is any approxi- 
mation whatsoever to A n , then the number 

xUc+l) = 7r n — tan -1 (yX^ k) ) (6-4) 

is a better approximation, as illustrated by Figure 6(b). The following 
theorem proves that this method always works. 


^or example, the buckling of a beam built in at one end and clamped at the 
other is governed by the equation 



tan 


/ Pi/ 

/ El ' 


where P is the axial buckling load, L is the length of the beam, E is Young's 
modulus and I is the usual moment of inertia. 
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Theorem: For arbitrary , the sequence { X^) converges to x n . 

The convergence is geometric with rate 

______ r 1 1 

Kn Ti<2n-1) yiT 2 (2n-l) 2 

Proof: For fixed y and n, let 

g(X) = nir - tan -1 (yX) 

and let I n denote the closed interval [(n-y)TT, (n+|-)TT]. Then gfX^ 0 )) e I n / 

regardless of the initial value X^°) . Hence all elements of the sequence 
{ X^ k ^ }y L= i belong to I n - It is sufficient to show that g is a contraction 
mapping on I n , because the iteration defined by (6-4) is just the usual 
Picard approximation scheme applied to the fixed point equation 


Xn _ 9 (X n ) . 

Since g is differentiable on I n , it suffices to show that 

sup { | g ' ( X ) | } < 1 . 

Xel„ 


But 


Therefore 


g r (X) 


yX 


l+y z X 


1+y 2 X^ X 


y 2 X 2 1 

l+y 2 X 2 yX 2 ’ 


and 


sup Ik'ixii) <_i xY P 1^) -Tiiir 
£F n ’ ,9 ' 1X> 11 - Xs?„ < TyX 7 T 1 ’ Y. 2 (2 4 „-l)^ 


QED 


The above algorithm has been coded as a function subroutine named 
PICARD using 


X ( Q) 


nTT 


as an initial value, and using the relative error test 

| X (k + ,) _ A <k)| < e | A (k+u| 


( 6 - 6 ) 


(6-7) 


for convergence . 


2 The relative error test (6-7) is useful whenever the exact limit is not 
known in advance and whenever the limits may be extremely large or extremely 
small. It does not guarantee that convergence is achieved to within e of 
the exact value, but is a reasonably general test and has worked well in 
the present study. 
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Typically, for n = 50, ten decimal convergence is achieved in 3-4 
iterations. The convergence characteristics of PICARD are shown in Table 1 
for several combinations of n and e in terms of the numbers 

k max( n ' e > = _<S a Y x <co {max{k+l: \^ +1) - ^ k) | < e |x^ k+1 )|}} (6-8) 

which represent the maximum number of iterations required for all values of y. 


Table 1. Convergence of function subroutine PICARD 
n = eigenvalue no. , k max = max no. iterations required 
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The rapid convergence with large n is explained by the convergence rate 
bounds given by equation (6-5) . Values of kmax were estimated by comparing 
results for 25 values of y which, to within the single precision accuracy 
of a CDC 6400, were given by values of y£ such that: 

i C.TT 

tan -1 Y£ = — ; -12 <_ £ <_ 12. 

Only one iteration is required for SL = 0 because the choice (6-6) of the 
initial value is then the exact eigenvalue, but the maximum number of itera- 
tions usually occurred for 1 = 1. This behavior can be understood from 
the graph of the eigenspectra vs. y. Since y can be any real number, the 
domain can be mapped onto a finite interval by employing the inverse tangent 
transformation. The resulting eigenspectra for all real values of y are 
shown in Figure 7. A gradually steepening shoulder with increasing n is 
observed near y = 0, which is where the maximum number of iterations should 
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be expected. This shoulder is merely exhibiting the fact that 


lim A n 


fa (n+y) 
■ Trn 

l it (n-— ) 


y < o' 
Y - 0 
y > o 


which is obvious from Figure 6(a). 

The infinite series (2-51) and its derivative are the result of a 
uniformly but slowly convergent series which has been modified for improved 
convergence. Consider the infinite series 


oo Otn 

f(6) = Z — exp(-A n 6), (6-10) 

n 1 ^n 

f ' (6) = - Z a n exp (-X n <5) , (6-11) 

n= 1 

where ot n and A n are given by (2-52) and (2-53) . From (6-9) , it follows 
that if y > 0 the terms of the series (6-10) and (6-11) for large n approach 
the corresponding terms of the series 

f(6) = Z exp(-Tr (n-^-) 6) , (6-12) 

n=i . 1. 2 

ir(n-j) 

f'(6) = - ? exp(-n(n-i-) 6) , (6-13) 

n=l 2 

respectively. The series (6-10) - (6-13) diverge if 6 <_ 0 and converge uni- 
formly on any closed interval such that 


6 > 0 . 


(6-14) 


The derived series for f' and f* are the more slowly convergent, converging 
geometrically with rate which approaches 1 as 6 approaches 0. However, 
the minimum argument needed for the collocation method is 


~ _ . |x-£ 1 

min 8n H ' 

where x is a collocation point and £ is a quadrature node 


(6-15) 


Although Hsu ' s 

interdigitation procedure tends to maximize this difference, the argument 
can still be quite small. Taking the number of quadrature points equal 
to the number of pressure basis functions. 


NP = NQ 


and referring to formulas (3-24) and (3-25) 

I 2i7T COS2 j TT | 7T 

l<i?3<NP ' COS 2NP+1 + 2NP+1 ' ~ C ° S 2NP+1 


one finds that 

2lt __ 3 tt 2 
COS 2NP+1 ~ 8NP 2 ’ 


(6-16) 


(6-17) 
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Consequent ly, applying a relative convergence criterion to the series 
(6-10) to (6-13) , it can be estimated that the number of terms NT required 
to achieve ND decimal convergence is given by 

NT « (NP) 2 - (6-18) 

For example if M = 0, Hh = 15, ND = 10 and NP = 10, then NT - 3000 terms. 

The series (6-10) and (6-11) appeared in the original expression for 
Bland's kernel. He improved their weak convergence properties by subtracting 
the corresponding terms of the series (6-12) and (6-13) and adding their 


closed form sums 

f (6) = ^ log coth -~- f (6-19) 

f ' (6) = -2 csch y- (6-20) 

to obtain 

f(<5) = Z [— - exp(-X n 6) - — exp (-7T (n-^-) 6) ] + — log coth (6-21) 

n_1 *n 77 (n-j) " * 4 

f ' (6) = [a n exp (-X n <5) - exp (-it (n-y) 6) ] - 2 csch yp (6-22) 


The infinite series appearing in (6-21) and (6-22) equal F and F 1 as defined 
by (2-51). F and F' are nonsingular at 6 = 0, whereas f, f * f f and f 1 are 
all singular at 6 = 0. The resulting smooth behavior of F and F' vs. 6 is 
depicted in Figure 8 for k = 0 and for y = 0^1 and 00 . In the case y = k = 0, 
it is possible to sum the series for F and F' in closed form: 

F(<5) = ~ " log ( 1+exp ( ) if y = k = 0, (6-23) 

F’ (6) = 7 if Y = k = 0, (6-24) 

. 7T O . 

1+exp (-—) 

thereby permitting an exact and independent check of numerical computations. 
Table 2 shows a comparison for small values of 6, including 6=0. 
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Table 

2 . Accuracy of 

subroutine SUM 

for y = k = 0 

(using a 

relative convergence test on F 1 

with e = 10” 10 ) 

6 

Decimals < 

of accuracy 

no. of terms 


F(5) 

F« (6) 


.OOO 3 

5 

0 

60057 

.001 

10 

6 

5496 

.002 

10 

7 

2859 

.003 

10 

8 

1485 

.020 

10 

9 

323 

.200 

10 

10 

37 


Thus, while the series reformulations described above have succeeded in 
removing the singularity and improving the convergence somewhat, it is 
apparent from Figure 8 and Table 2 that additional gains in efficiency 
remain to be accomplished. Variations in frequency and Mach number are 
shown in Figures 9-11. 

The prospects for improving the efficiency of subroutine SUM appear 
favorable for several reasons. The present version is accurate and pro- 
vides a sound basis for comparison. The data presented above represent a 
worst case in the sense that convergence is more rapid with increasing y 
(y = 0 corresponds to an open tunnel, y = 00 to a closed tunnel) . 


3 In the case y = k = 0, the series for F(0) is an alternating harmonic 
series and converges slowly to 

F (0) = log 2, 

but the derived series is null and does not converge to 

F' (0) = j. 

This is due to the non uniform convergence of F(5) at 6 = 0 so that inter- 
change of summation and differentiation is not valid. 
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§7. Use of the TWODI program 

This section describes the preparation of input data for the pilot 
version of the TWODI program. It is assumed that the user knows the physical 
meaning of all input and output quantities, and the procedure for submitting 
a run at the computer facility. We reiterate that the solution is based 
upon the assumption of linearized inviscid subsonic potential flow about 
a thin airfoil located midway between two parallel walls using the boundary 
condition (2-49) . 

TWODI will predict any combination of pressure, section coefficients 
and generalized aerodynamic forces for the following primary parameters: 

( 1 ) mode shape , 

(2) reduced frequency, 

(3) Mach number, 

(4) tunnel depth to airfoil chord ratio, 

(5) tunnel wall ventilation coefficient. 

Additional parameters are the number of pressure terms used in the solution 
process, the number of points at which pressure is to be calculated, etc. 

TWODI will operate in either BATCH (noninteractive card jobs) or TIMESHARE 
(interactive remote terminal) modes. Input and output are fully compatible 
between BATCH and TIMESHARE. A problem is defined by a set of data sufficient 
to cause execution of the program, and contains from 1 to 50 combinations 
of (M, k, r^, c^) , called data cases . As many problems may be loaded in succession 
as desired. All input data and their acceptable values are precisely defined 
in Section 7.1 below. Input data are automatically checked for acceptabi- 
lity and unacceptable data cases or unacceptable problems will be select- 
ively deleted by TWODI with an explanatory comment. The input format is 
defined in Section 7.2 consisting of six input data units. For BATCH runs 
the format may be formal on a column by column basis, or it may be free 
with individual data separated by a comma or blanks. Each input data unit 
(and certain subunits) must begin on a new card. TIMESHARE input is prompted 
in a self-explanatory fashion, and the interactive communication between 
TWODI and user is described in Section 7.3. At the present time, compre- 
hensive estimates of computer time and the number NP of pressure modes 
required to achieve a certain accuracy are not available. 




7.1 Definition of input data 

TITLE - From 1 to 72 alphanumeric characters of user-selected title. 

Printed at the top of each page, with the time and date. 

NC - Number of M,k,riH' c W data cases to follow for this problem. 

Each problem consists of NC cases. As many problems may be 
run as desired. 1<NC<50 

NP - Number of aerodynamic pressure modes to be used in the solution 

process. 1<NP<30 

NH - Number of airfoil deflection modes. 1<NH<5 

NX - Number of points at which the airfoil deflection modes are to 

be collocated to the input data. 1<NX<10 

NL - Number of loading points at which pressure is to be calculated. 

°<NL<5° 

SEC - Logical variable whose truth value is to calculate section 

coefficients (lift, pitching moment, center of pressure) for 
each deflection mode. T or F 

GAF - Logical variable whose truth value is to calculate the general- 

ized aerodynamic force coefficient matrix. T or F 

TSH - Logical variable whose truth value is interactive timesharing 

remote terminal operation. T or F 

XM(M) - Chordwise coordinate of matching point M of NX for airfoil 

centerline deflections. Nondimensionalized by the semichord, 

-1 at the leading edge, +1 at the trailing edge. -1<_XM(M)<1. 

If M / N, XM (M) ± XM(N) . 

HM(NR,M) - Vertical coordinate of matching point M for airfoil centerline 
deflection mode NR. Nondimensionalized by the semichord, posi- 
tive up. 

XL (L) - Chordwise coordinate of loading point L of NL at which pressure 

is to be calculated. To be omitted if NL = 0. -1<XL(L)<+1 

MACH - Mach number, M. C<_M<1 

FREQ - Airfoil reduced frequency, k, referred to semichord. 

CW - Tunnel wall porosity coefficient, cw- (Xc^ 00 . Use c w = 10 100 

for closed wall conditions. 

ETAH - Tunnel semiheight nondimensionalized by the airfoil semichord. 

Equals the tunnel height to airfoil chord ratio nH- Positive. 
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7.2. Input data format 

DIGITAL COMPUTER INPUT DATA SHEETS 
Programmer pate Page of Ident . 


DATA 

SEQUENCE 

DESCRIPTION - DO NOT KEYPUNCH 

l 

TITLE 

INPUT 
DATA UNIT 

I 

7 3 SO 

VARIABLE TITLE 

1 3 


2 5 

72 ALPHANUMERIC CHARACTERS 

3 7 


4 9 


6 1 


1 NC NP 

INPUT 
DATA UNIT 

II 

73 SO 

CONTROL PARAMETERS 

1 3 

NH NX 


2 5 

NL SRF 

SEPARATE INDIVIDUAL DATA 

3 7 

GAF TSH 

WITH A COMMA OR BLANKS 

4 9 


6 1 


1 XM (1) 

INPUT 
DATA UNIT 

III 

7 3 6 0 

MATCHING POINTS 

1 ! 

XM(2) 


2 5 

SEPARATE INDIVIDUAL DATA 

3 7 

XM (NX) 

WITH A COMMA OR BLANKS 

4 9 


6 1 


1 HM (1, 1) 

INPUT 
DATA UNIT 

IV 

7 3 6 0 

MODE SHAPES 

HM(2 , 1) 


2 5 

SEPARATE INDIVIDUAL DATA 

J 7 

HM (NX, 1) 

WITH A COMMA OR BLANKS. START 

4 9 

EACH MODE ON A NEW CARD 

6J 





EACH MODE ON A NEW CARD 
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Programmer 


DIGITAL COMPUTER INPUT DATA SHEETS 
Date Page of Ident . _ 


SEQUENCE 


DESCRIPTION - DO NOT KEYPUNCH 


XL(1) 


XL (2) 


XL(NL) 


INPUT 
DATA UNIT 


LOADING POINTS 


OMIT IF NL = 0 


SEPARATE INDIVIDUAL DATA 


WITH A COMMA OR BLANKS 


1 

MACH 

1 3 

FREQ 

2 5 

ETAH 

3 7 

CW 


INPUT 

FLOW PARAMETERS 

DATA UNIT 


VI 

SEPARATE INDIVIDUAL DATA 


WITH A COMMA OR BLANKS. START 


EACH CASE ON A NEW CARD 

7 3 8 0 



REPEAT INPUT DATA UNITS I-VI 


FOR AS MANY ADDITIONAL 


PROBLEMS AS DESIRED 
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7.3. Prompted interactive input 

The following is an example of automatically prompted interactive 
communication between TWODI and its user. Responses by the user are 
underlined. The particular data are for the problem of Section 10, and 
the output is presented as the APPENDIX. 


ENTER TITLE, 

? C0MPAR I SON WITH THE KUESSNER -SCHW ARZ SOLUTION USING ETAH = 300 
ENTER NC,NP,NH,NX,NL,SEC,GAF,TSH 
? l,10 y 5,5,20 > T.T.T 
ENTER MATCHING POINTS, 

5,0, .5,1 

ENTER MODE SHAPE FOR MODE 1 

? 1 , 1 , 1 , 1,1 

ENTER MODE SHAPE FOR MODE 2 

9_ 3 f _2 , -1 , 0 , 1 

’enter mode shape FOR MODE 3 
? 5, 1,-1, -1,1 

ENTER MODE SHAPE FOR MODE 4 

? - 7 , 1 , 1 , - 1,1 

ENTER MODE SHAPE FOR MODE 5 
? 9, -2, 1,0,1 

ENTER LOADING POINTS 

9, -.8, -.7, -.6, -.5, -.4, -.3, -.2, -.1,0, •l,«2,.3,.4,.5/«6,.7/.8 

ENTER MACH, F R EQ,ETA H , CVrFO R^CAS E~1 " 

?0,1,300,1E100 


• 9,1 
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§8. Convergence characteristics of TWODI 


Numerical calculations with the TWODI program reflect the theoretical 
convergence predicted mathematically in Section 4. Tables 3 and 4 show 
for M = 0 and k = 0 the convergence of lift coefficient and center of pressure 
vs. NP with orders of magnitude variation in tunnel depth to chord ratio 
and for both the open and closed wall conditions. For the open wall condi- 
tion, slower convergence of the infinite series for F and F* resulted in 
our omitting some of the calculations for the deeper tunnel cases. 


Table 3. Convergence of CL a vs. NP for M = 0 and k = 0 


(flat plate at unit angle of attack) 


Wall 

NP 

a 

n 

n H =l0 

riH=100 

riH=1000 

n H =ioooo 

Open 

1 

1.91357 

5.39195 

6.18551 

6.27333 

6.28220 

H 

2 

1.91357 

5.39195 

6.18551 

6.27333 

6.28220 

ft 

3 

1.91357 

5.39195 

6.18551 

6.27333 


It 

4 

1.91357 

5.39195 

6.18551 



II 

5 

1.91357 

5.39195 




It 

10 

1.91357 





Closed 

1 

9.20520 

6.30906 

6.28344 

6.28319 

6.28319 

ii 

2 

8.26642 

6.30894 

6.28344 

6.28319 

6.28319 

it 

3 

8.30090 

6.30894 

6.28344 

6.28319 

6.28319 

it 

4 

8.29957 

6.30894 

6.28344 

6.28319 

6.28319 

it 

5 

8.29957 

6.30894 

6.28344 

6.28319 

6.28319 

n 

10 

8.29957 

6. 30894 

6.28344 

6.28319 

6.28319 


Table 4. 

Convergence of x^p 

vs. NP for 

M = 0 and 

o 

n 


(flat plate 

at unit 

angle of attack) 


wall 

NP 

n H =i 

n H =l0 

riH=100 

tih=1000 

riH=i0000 

Open 

1 

.250000 

.250000 

.250000 

.250000 

.250000 

ii 

2 

.110850 

.247953 

.249979 

.250000 

.250000 

i» 

3 

.111373 

.247954 

.249979 

.250000 


ii 

4 

.111440 

.247954 

.249979 



m 

5 

.111434 

.247954 




it 

10 

.111435 





Closed 

1 

.250000 

.250000 

.250000 

.250000 

.250000 

n 

2 

.307484 

.251020 

.250010 

.250000 

.250000 

ii 

3 

. 306176 

.251019 

.250010 

.250000 

.250000 

tt 

4 

.306172 

.251019 

.250010 

.250000 

.250000 

n 

5 

.306175 

.251019 

.250010 

.250000 

.250000 

ii 

10 

.306175 

.251019 

.250010 
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Whereas section coefficients depend upon only the first two Fourier 

coefficients P^ and P^, the generalized aerodynamic forces and the details 

of the pressure functions depend upon all the Fourier coefficients P,n ; 

r n 

n = 1,...,NP. Tables 5 and 6 show, for steady and unsteady flow respectively, 
the convergence of these truncated sequences of Fourier coefficients with 
respect to the number NP of pressure basis functions employed, for both the 
open and closed tunnel wall conditions. The convergence for steady flow 
is much more rapid than for unsteady flow. 


Table 5. Convergence of P^ vs. NP for M = 0, k = 0 and Hjj = 



(flat plate 

at unit angle 

of attack) 


Wall 

NP 

p ^i 


P ^3 

P ^4 

Open 

1 

3.43262 




ii 

2 

3.43262 

-.014052 



■I 

3 

3.43262 

-.014047 

.000023 


ii 

4 

3.43262 

-.014047 

.000023 

.000006 

n 

5 

3.43262 

-.014047 

.000023 

.000006 

Closed 

1 

4.01647 




ti 

2 

4.01640 

.008194 



it 

3 

4.01640 

.008188 

-.000023 


ii 

4 

4.01640 

.008188 

-.000023 

-.000006 

ii 

5 

4.01640 

.008188 

-.000023 

-.000006 


Table 

6 . 

Convergence of P^ n 

vs. NP for M = .5, k = .1 and r) H = 10 

(flat plate 

oscillating about 

midchord-unit maximum angle of attack) 

Wall 

NP 

p >h 

P 4 , 2 P 4»3 

Open 

1 

3 . 603 08- . 26131 Oi 


it 

2 

3.64729-. 2603 06i 

-.007670+. 522590i 

ti 

3 

3. 64763- . 260345i 

-.007831+. 523485i -. 010920+. 0000241 

ii 

4 

3 . 64772- . 260359i 

-.007871+. 523535i -. 011029+. 000015i 

M 

5 

3 . 64776- . 260365i 

-.007853+. 523553i -.011068+. 000013i 

ii 

10 

3. 64779-. 2603711 

-.007898+. 523572i -. 011106+. 000012i 

H 

15 

3. 64780- . 2 60372 i 

-.007900+. 523574i - . 011111+. 000012i 

Closed 

1 

3.77355-. 7622 77 i 


n 

2 

3 . 81694- . 772065i 

.038394+. 519791i 

tt 

3 

3 . 8172 9-. 7722 03i 

. 038279+ . 520718i -.010978+. 0008551 

ii 

4 

3. 81739- . 7 72 245 i 

. 038244+ . 520744i - . 011093+. 0008611 

n 

5 

3 . 81743- . 77 2 2 62 i 

. 038232+. 520794i -.011134+. 000865i 

it 

10 

3 . 81747- . 772280i 

. 038220+. 520815i -. 011174+. 000869i 

n 

15 

3.81747-.772282i 

. 038219+ . 5208l7i 011787+. 000869i 



In general, we have found that the convergence of the TWODI program 
with respect to the number of pressure basis functions is remarkably good 
for steady flow and good for unsteady flow with the program being most 
efficient for narrow tunnels and low frequencies- The practical utility 
of the TWODI program would be improved by the incorporation of more 
efficient computer algorithms for deep tunnels and higher frequencies. 
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§9. Comparison with the Sohngen solution 


For steady incompressible flow in an infinite atmosphere, the results 
of the TWODI program can be compared against the exact closed form Sohngen 
solution. This solution is given by the Sohngen inversion formula 1 (2-25) 
or, equivalently, by Bland* s solution (3-24) in terms of airfoil polynomials. 
We shall utilize the latter because of their elegant simplicity. 

Since the downwash polynomials are complete, they may be used as basis 
functions for deflections. In the present comparison, airloads will be 
calculated for airfoil deflections, or contours, spanned by the first five 


downwash polynomials (Fig. 3) : 

h n = Xn' n = 1,2, 3,4,5. (9-1) 

Written out, the first five downwash and pressure polynomials are 

Xl (x) = 1, (9-2 i ) 

X2 (*) = ~l+2x, (9-22 ) 

X 3 (x) = -l-2x+4x 2 , ( 9—2 3 ) 

X4 (x) = 1-4 x-4x 2 +8x 3 , ( 9-2 4 ) 

X5<x) = l+4x-12x 2 -8x 3 +16x 4 , (9-25) 

<h (?) = 1, (9-3 1) 

iM£> = 1+25, (9-3 2 ) 

<|*3 (C) = -1+2C+4? 2 , (9-3 3) 

<1*4 (5) = -l-45+45 2 +85 3 , (9-3 4) 

<J* 5 (?) = l-45-125 2 +85 3 +165 4 . (9-3 5 ) 

Then the downwash functions 

w n = xA (9-4) 


Extensive closed form integrals based on the Sohngen and Kussner-Schwarz 
solutions have been tabulated by Fromme [88, 1964,App. B and C] for piecewise 
continuous downwash functions expressed as powers. 
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are found to be 

wj = 0, (9-5!) 

w 2 = 2 * Xl f (9-5 2 ) 

W 3 = 2 xi+4x 2 / (9-5 3 ) 

w 4 = 4xi+2x 2 +6X3» (9-5 4 ) 

w 5 = 4 xi+6x2+ 2 X3 + 8 X4/ (9-5 5 ) 

and by equation (3-24 ) , the corresponding pressures are 

Api (5) = 0 , (9-6!) 

Ap 2 (S) = - / IT§ [8+1 <5) (9-6 2 ) 

Ap 3 (£) = " [8*i(e)+16* 2 (5)], O- 63 ) 

Ap 4 (C) =- [16^! (5)+8i|- 2 (C)+24i(»3(C)], (9-6 4 ) 

Aps (£) =- /l^| t 16 ’!'! (0+24^2 (£)+8^3(C)+32 i|/ 4 (?)]. (9-6 5 ) 


The results in (9-6) , based on the Bland transform, have been checked in- 
dependently with the formulas of Fromme, and are in agreement. Numerical 
values of pressure are presented in Table 7. Figure 12 depicts the deflection 
basis functions and the corresponding pressure functions. These results 
have been used to check the predictions of the TWODI code using six pressure 
basis functions with nn - 10 4 and c w = 10 100 . To within the six decimal 
accuracy of printout, 2 the pressures predicted by TWODI are correct for 
all five deflections. 


2 In this work, we adopt a standard of six decimal accuracy as a compromise 
between reasonable generality and practical utility. Greater accuracy is 

attainable, but we note that six decimal accuracy is already orders of 
magnitude above experimental accuracy. However, intermediate calculations, 

especially those which occur many times, are performed to higher accuracies 
that are established internally in the computer program. 
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Table 7. Pressures for the Sohngen comparison 



A Pl ( 5 ) 

AP 2 (?) 

Ap 3 (?) 

Apit (£) 

Aps (?) 

-.9 

0.0000 

- 34.8712 

20.9227 

- 87.8754 

- 2.51073 

-.8 

0.0000 

- 24.0000 

4.80000 

- 30.7200 

- 67.5840 

-.7 

0.0000 

- 19.0438 

- 3.80876 

- 5.33227 

- 84.2498 

-.6 

0.0000 

- 16.0000 

- 9.60000 

7.68000 

- 81.4080 

-.5 

0.0000 

- 13.8564 

- 13.8564 

13.8564 

- 69.2820 

-.4 

0.0000 

- 12.2202 

- 17.1083 

15.6419 

- 53.1823 

-.3 

0.0000 

- 10.9022 

- 19.6239 

14.3909 

- 36.3696 

-.2 

0.0000 

- 9.79796 

- 21.5555 

10.9737 

- 21.0068 

-.1 

0.0000 

- 8.84433 

- 22.9953 

6.01415 

- 8.56131 

0 

0.0000 

- 8.00000 

- 24.0000 

.000000 

.000000 

.1 

0.0000 

- 7.23627 

- 24.6033 

- 6.65737 

4.11020 

.2 

0.0000 

- 6.53197 

- 24.8215 

- 13.5865 

3.55339 

.3 

0.0000 

- 5.87040 

- 24.6557 

- 20.4290 

- 1.54978 

.4 

0.0000 

- 5.23723 

- 24.0913 

- 26.8146 

- 10.7258 

.5 

0.0000 

- 4.61880 

- 23.0940 

- 32.3316 

- 23.0940 

.6 

0.0000 

- 4.00000 

- 21.6000 

- 36.4800 

- 37.2480 

.7 

0.0000 

- 3.36067 

- 19.4919 

- 38.5805 

- 51.0016 

.8 

0.0000 

- 2.66667 

- 16.5333 

- 37.5467 

- 60.7573 

.9 

0.0000 

- 1.83533 

- 12.1132 

- 31.0537 

- 59.3324 

1.0 

0.0000 

0.00000 

0.00000 

0.00000 

0.00000 


Values of section coefficients may be obtained by inspection of (9-6) 
using (3-37) and (3-38) . The results are given in Table 8. Again the 
predictions of TWODI are correct to within the six decimals of printout. 


Table 8. Section coefficients for the Sohngen comparison 


Mode 

CL 

Cm 

1 

0 

0 

2 

-4tt 

0 

3 

-4tt 

-4tt 

4 

-8*jt 

-2t T 

5 

-8tt 

-6tt 


To calculate the generalized aerodynamic force matrix, we expand the 
deflection functions given by (9-1) in terms of the pressure polynomials. 
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It is easy to show that 

XI = ♦if 
X 2 = -^l+^Z , 

X3 = 2if) 1 -2^2+<|'3, 

X4 = -2 ’l'l + 2’f , 2~ 2 'i J 3 + 'J'4r 

X 5 = i — 2ip2 + 2i/J 3 — 2 ^ 4 +ip 5 . 

Combining equations (9-6) , (9-7) and (3-42) , we have our final result. 


0 

-4tt 

-4tt 

-8tt 

-8tt 

0 

8tt 

0 

12tt 

4tt 

0 

-8tt 

8tt 

-20tt 

4tt 

0 

8tt 

-8tt 

32tt 

-16tt 

0 

-8tt 

8ir 

-32tt 

32tt 


Again, the TWODI program is correct to within the accuracy of printout. 


(9-7 j ) 
(9-7 2 ) 
(9-7 3 ) 
(9-7 4 ) 
(9-7 5 ) 


(9-8) 




§10. Comparison with the Kussner-Schwarz solution 


n 

This section compares the exact closed form Kussner-Schwarz solution 
for unsteady incompressible flow in an infinite atmosphere against the 
results of the TWODI program. Also, this section presents new and simpli- 
fied expressions using airfoil polynomials for unsteady airloads based 
on the Kussner-Schwarz solution. 

For airfoil deflections spanned by the first five downwash polynomials, 
the downwashes are 

w n = + ik) Xn » n = 1,2, 3, 4, 5 (10-1) 

and it follows that 


wi = ikXlf (10-2J) 
w 2 = 2xi+ikx 2 , (10-2 2 ) 
w 3 = 2 Xl+4x 2 + ikx 3 , (10-2 3 ) 
w 4 = 4xi+2x 2 +6x3+ikx 4 , (10-2 4 ) 
w 5 = 4xi+6x 2 +2x 3 +8x4+ikx5. (10-2 5 ) 


To organize the calculations, we make use of the inverse of the Bland trans- 
form, the Lambda transform and the Theodorsen functional, defined as follows: 


H-lw(C) = i / 1 /^ 

7T _ 1 v 1 — X X-C 


l t 1 

Aw (£) = — y I A (x, £ ) w (x) dx 


Hw = - / 


/ 1+x 


IT -i / l-> 


w(x) dx 


Thus, 


A p(£) =4 y - ikA + [l-c(k)]f/}w(5). 


(10-3) 

(10-4) 

(10-5) 

(10-6) 


(a) The inverse Bland transform. Since 
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then 

H _1 Wx = ik^i/ (10-8!) 

H -1 W2 = 21^+ikiJj! , (IO-82) 

H“ 1 W3 = 2i(;i+4\jj2+ik^3 , (IO-83) 

H _1 W4 = 4i/; 1 +2^2+6^ 3 iki(y4 , (IO-84) 

H -1 W5 = 4i|; 1 +6^2+ 2l i J 3 + 8^4+ik^5. (IO-85) 


(b) The Lambda transform , Fromme 1 has shown that 


Hence, 


A [1] 

= 1 + Cf 


(10-9!) 

A[x] 

- ^f 52 - 


(10-9 2 ) 

A [x 2 ] 

■ ? + I 5 + I 52 + 

if 3 
3T ' 

(10-9 3 ) 

A [x 3 ] 

= h + h 2 + V 

+ V 

4* ' 

(IO- 94 ) 

A[x 4 ] 

3 J. -3_r . JL, 
= 40 + 40* + 10 5 

2 + _i_f3 4. if 4 . if 5 

+ 10 4 I 4 + s 4 ' 

(10-9 5 ) 

Axi = 

^1 + 


(10-10!) 

AX2 = 

1, 1, 1, 
f’/'l “ ^2 + ^3' 


(10-10,) 

Ax 3 = 

1 , 1 , 

- 4*2 - Y^3 

1 . 

+ ^ 4 , 

(IO-IO 3 ) 

Ax 4 = 

" 6 *’ " 

1 , 1, 

+ 8^5 ' 

(IO-IO 4 ) 

Ax 5 = 


i ^ 4 ~ ^ + Y^ 6 - 

(10-10 5 ) 


1 [ 88 , 1964 , Appendix C, formulas (C-45) to (C-49) ] . 
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Combining equations (10-2) , (10-3) and (10-10) , we get 


Awj = jiki 

-ik’t/2 f 

( 10 -llx) 

Aw 2 = (l-“ik)ipi 

+ (l--|ik )^ 2 + jik^ 3 . 

( 10 - 11 2 ) 

AW 3 = - i|j i 

- ^ ik 'l '2 + (1-J^ik)ij>3 + 

(IO-II 3 ) 

Aw 4 = \p i 

- ^ik'1'3 + (l-^-ik)^4 + 

(IO-II 4 ) 

AW 5 = - ^ 1 

- |ik ^4 + (l-^lk )* 5 + ^±k* 6 . 

(IO-II 5 ) 

(c) The Theodorsen functional. Since 



n {|/l2 Xm(x)x n (x)dx = 6 ron 


and since Xl (x) = 

l r it follows from (10-2) and (10-5) that 



H W]_ = ik. 

( 10 - 12 !) 


H W 2 = 2 , 

( 10 - 12 2 ) 


Ww 3 = 2 , 

(IO-I 23 ) 


Hw 4 = 4 , 

(IO-I 24 ) 


ff W 5 = 4. 

( 10 - 12 5 ) 

The Fourier 

solution for pressure now follows. 


Api (?) 

= /^| ( [-4ikC(k)+2k 2 ]0i (?)+2k 2 ^ 2 (?)>. 

(10-13!) 

Ap 2 (?) 

= / ( t- 8 C (k) -4ik-2k 2 ] i> l (?) - (k 2 + 8 ik) 1> 2 (?) +k 2 t p 3 (?) }, 

(10-13 2 ) 

Ap 3 (?) 

= / { [~ 8 C (k) +4ik] i/ii (? ) - (16+k 2 ) i /) 2 (? ) 



-( 8 ik+Yk 2 )* 3 (?)+|k 2 K(?)}/ 

(10-13 3 ) 

Ap 4 (?) 

= /’j^ { [-16C(k)-4ik]*i (?) - 8^2 (?) “ (24+|ic 2 )^3 (?) 



- ( 8 ik+|k 2 ) *4 ( ? ) +|k 2 i/j 5 (? ) , 

(10-13 4 ) 

A P 5 (?) 

= / { [-16C (k) +4ik] !}<! (?) —24^2 (?) ~8^ 3 (?) 



- (32+|k 2 ) i|<4 (? ) - (8ik+^k 2 ) <P 5 (? ) +|k 2 i|/ 6 (? ) } . 

(10-13 5 ) 
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The section coefficients follow by inspection of (10-13) 

C L = IT [-2ikC(k)+k 2 ] , (10-14]) 

C L = iT[-4C(k)-2ik-k 2 ] , (10-14 2 ) 

2 

C L = 7r[-4C(k)+2ik] , (10-14 3 ) 

3 

Cl = TT[-8C(k)-2ik] , ( 10-14 4 ) 

4 

Cl^ = tt [-8C (k) +2ik] , (10-14 5 ) 

C Ml = TT(|k 2 ) , (10-15]) 

c M 2 = 17 (^k 2 -2ik) , (10-15 2 ) 

C M = tt(— 4— ~k 2 ) f ( 1 0—15 3 ) 

= it (-2) , (10-15 4 ) 

Cm 5 = T! (-6) , (10-15 5 ) 


and the generalized aerodynamic forces follow from (3-42) using (9-7) and 
( 10 - 11 ) . 
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An = 

7 r[- 2 ikC(k)+k 2 ], 

(10-16!) 

a 12 = 

u[-4C(k)-2ik-k 2 ] , 

(10-16 2 ) 

a 13 = 

it [-4C (k) +2ik] , 

(IO-I 63 ) 

a 14 = 

tt [- 8 C (k) -2ik], 

(IO-I 64 ) 

Al5 " 

it [- 8 C (k) +2ik] , 

(IO-I 65 ) 

a 21 * 

ir[4ikC(k)-k 2 ] , 

(10-16g) 

A 2 2 = 

tt l 8 C(k)+|k 2 ] , 

(IO-I 67 ) 

a 23 = 

Tr[8C(k)-8-4ik-|k 2 ], 

(10-16g) 

A 2 i+ “ 

tt [16C (k) -4+4ik] , 

(10-16g) 

a 25 - 

tt [16C (k) -12~4ik] , 

(10-16 ! 0 ) 

A 3l » 

77 [-4ikC (k) 3 , 

(10-16 j ! ) 

A 3 2 = 

7 t [-8C (k)+4ik-|k 2 ] , 

(10-16 12 ) 

A 3 3 = 

tt [- 8 C (k) +16+^k 2 ] , 

(10-16 23 ) 

A 3 4 = 

n [-16C(k)-4-4ik-|k 2 ], 

(10-16! 4 ) 

A 35 » 

Tr[-16C(k)+20+4ik] , 

( 10 - 16 ! 5 ) 

A 41 = 

tt [4ikC (k) ] , 

(10-16i 6 ) 

II 

CM 

4 

tt [ 8 C (k) - 4ik] , 

(10-16 17 ) 

A 43 = 

tt [ 8 C (k) -16+4ik-|k 2 ] , 

(10-16 18 ) 

A 44 = 

7r[16C(k)+16+^k 2 ], 

(IO-I 61 9 ) 

A45 = 

ir[16C(k)-32-4ik~|k 2 ] , 

(10-16 20 ) 

a 5 i = 

it [-4ikC (k) ] , 

(10-16 2 i ) 

a 52 = 

it [-8C(k)+4ik] , 

(10-16 22 ) 

A53 = 

tt [- 8 C (k) +16-4ik] , 

(10-16 23 ) 

A54 = 

tt [- 16C(k)-16-4ik-|k 2 ] , 

(10-16 z4 ) 

A55 = 

it [-16C(k)+48+^k 2 ] . 

(10-16 25 ) 
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Numerical values of airloads based on the exact values given by (10-13) 
to (10-16) are presented in Tables 9 to 11 and Figure 13. The pilot version 
of TWODI will not accept Oh = °°f and in its present form is most efficient 
for small values of riH and Thus, the Kussner-Schwarz comparison re- 

presents a worthy test. Numerical values predicted by TWODI for nn = 300, 
k = 1 and NP = 8, presented as the APPENDIX, are seen to be in close 
agreement. The aeroelastic effect is bounded by the relative error supremum 
norm of the generalized aerodynamic force matrix: 


e [A] = 


|Ag| ACT (n H =”) ~ Ags ODI (n H =300) H sup 


EXACT 


= .0859% 


■rs (Hh- 00 ) llsup 


(10-17) 


This error norm represents the maximum difference between the exact solution 
in an infinite atmosphere and the numerical solution by TWODI, considering 
the combined effects of numerical inaccuracies, finite tunnel depth, and of 
using only eight pressure basis functions. 
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Table 9. Pressures for the Kussner-Schwarz comparison with k - 1 



Re(Apj) 

Re(Ap 2 ) • 

Re(Ap 3 ) 

Re(Ap 4 ) 

Re (Ap 5 ) 

-.9 

-0.00476 

-22.1235 

39.8542 

-58.0293 

31.2030 

-.8 

1.19673 

-17.2664 

19.0216 

-10.3985 

-45.0000 

-.7 

1.90178 

-15.1291 

7.87586 

11.3355 

-67.3412 

- . 6 

2.39782 

-13.7510 

0.15838 

22.4864 

-68.0052 

-.5 

2.76939 

-12.6708 

-5.74258 

27.4860 

-58.0773 

-.4 

3.05338 

-11.7245 

-10.4536 

28.3354 

-43.2973 

-.3 

3.26916 

-10.8415 

-14.2720 

26.1859 

-27.2281 

-.2 

3.42795 

-9.98838 

-17.3565 

21.8174 

-12.2751 

-.1 

' 3.53653 

-9.14889 

-19.7975 

15.8246 

-0.10282 

0 

3.59891 

-8.31548 

-21.6488 

8.70238 

8.16904 

.1 

3.61715 

-7.48545 

-22.9421 

0.89173 

11.8675 

.2 

3.59169 

-6.65892 

-23.6943 

-7.19346 

10.7190 

.3 

3.52144 

-5.83773 

-23.9107 

-15.1462 

4.83414 

.4 

3.40348 

-5.02478 

-23.5855 

-22.5477 

-5.27874 

.5 

3.23253 

-4.22359 

-22.6988 

-28.9431 

-18.6663 

.6 

2.99945 

-3.43774 

-21.2084 

-33.8024 

-33.8236 

.7 

2.68808 

-2.66984 

-19.0296 

-36.4390 

-48.4592 

.8 

2.26630 

-1.91849 

-15.9772 

-35.8007 

-58.8993 

.9 

1.65154 

-1.16439 

-11.5352 

-29.6921 

-57.9949 

1.0 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 


c. 

Im (Apj ) 

Im (Ap 2 ) 

Im ( AP 3 ) 

Im (Ap 4 ) 

Im (Ap 5 ) 

-.9 

-9.40537 

13.9580 

5.58891 

-10.7213 

40.2743 

-.8 

-6.47322 

4.80655 

15.3666 

-23.1229 

41.3507 

-.7 

-5.13645 

0.00520 

19.8108 

-25.0513 

32.0497 

-.6 

-4.31548 

-3.19563 

21.7644 

-22.5833 

20.3991 

-.5 

-3.73731 

-5.53878 

22.1740 

-18.0058 

9.70705 

-.4 

-3.29600 

-7.32879 

21.5109 

-12.5557 

1.50242 

-.3 

-2.94050 

-8.71875 

20.0630 

-7.01504 

-3.63101 

-.2 

-2.64268 

-9.79529 

18.0309 

-1.91505 

-5.69315 

-.1 

-2.38547 

-10.6108 

15.5684 

2.37511 

-5.06828 

0 

-2.15774 

-11.1978 

12.8022 

5.60437 

-2.39563 

.1 

-1.95175 

-11.5761 

9.84331 

7.61651 

1.52646 

.2 

-1.76179 

-11.7558 

6.79503 

8.33836 

5.81963 

.3 

-1.58335 

-11.7392 

3.75865 

7.77561 

9.59778 

.4 

-1.41257 

-11.5205 

0.83939 

6.01520 

12.0401 

.5 

-1.24577 

-11.0839 

-1.84626 

3.23568 

12.4733 

.6 

-1.07887 

-10.3989 

-4.15891 

-0.26982 

10.4758 

.7 

-0.90643 

-9.40897 

-5.91387 

-4.04442 

6.03222 

.8 

-0.71925 

-7.99927 

-6.82594 

-7.34788 

-0.18414 

.9 

-0.49502 

-5.87254 

-6.31302 

-8.75717 

-6.07319 

1.0 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 
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Table 10. Section coefficients for the Kussner-Schwarz 
comparison with k = 1 


Mode 


c L 


C M 

1 

2.51156 

- 

3 . 38937i 

1.57080 

2 

-9.92033 

- 

5.023121 

-.785398 - 6 . 28319i 

3 

-6.77874 

+ 

7 . 54325i 

-13.3518 

4 

-13.5575 

- 

3 . 76305i 

-6.28319 

5 

-13.5575 

+ 

8. 80332i 

-18.8496 


Table 11. Generalized aerodynamic forces for the 
Kussner-Schwarz comparison with k = 1 


s 

H 

II 

U 

r - 2 

r = 3 

r = 4 

r = 5 

1 

2.51156 
-3 . 38937i 

-9.92033 
-5 . 02312i 

-6.77874 

+7.543251 

-13.5575 
-3 . 76305i 

-13.5575 

+8.803321 

2 

-1.88153 
+6 . 77874i 

18.2699 
-2 . 52013i 

-13.1461 

-15.08651 

14.5486 

+7.5261H 

-10.5842 
-17 . 6066i 

3 

-1.26007 

-6.778741 

-15.1283 
+15 . 0865i 

39.3260 
+2. 52013i 

-40.7285 
-7. 52611i 

35.7169 
+17 . 6066i 

4 

1.26007 

+6.77874i 

13.5575 
-15 . 0865i 

-37.7520 
+10 . 0462i 

79.2130 
-5 . 04027i 

-74.2014 

-17.60661 

5 

-1.26007 
-6 . 77874i 

-13.5575 
+15 . 0865i 

36.7080 
-10. 0462i 

-78.1658 
+17 . 6066i 

125.095 
+5 . 04027i 
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§11. Verification and extension of Bland 's results 


Eland [5 , 1970] presented numerical results for a flat plate oscillating 
about the 42.5% chord at M = .85 in a closed wind tunnel with dh = 7.5. 

His published values of lift coefficient magnitude for various frequencies 


ClI 

- 12.2351 

k = 

0 

c lI 

= 7.99420 

k = 

.1 

c l I 

- 5.43549 

k - 

.2 


have been verified with the TWODI program, as have his values for generalized 
forces, etc. 

Bland’s results were based on a computer program written only for the 
closed tunnel condition cyj = We have extended these results to include 

the effect of ventilating the wind tunnel walls. Table 12 shows the effect 
on Cl^ of variations in the parameters riH and cw and indicates continuous 
behavior with respect to both parameters. Values of c w 10 6 reproduce 

Bland’s closed wall results for nn = 7.5 and the lift coefficient for 
M = .85, njj - 1000 and c w = 10 6 agrees to six decimals with the well known 
infinite atmosphere solution C = 2 t t/£ 3. We have since observed that the 
numerics are not adversely sensitive to large values of cw and that the 
convergence rate of the infinite series in the kernel improves somewhat 
with large values of c^/hh- Consequently we now use c w = 10 100 to represent 
a closed wall condition. 

Table 12 also indicates that the effect of ventilation is considerably 
more pronounced for narrow tunnels than deep tunnels. Furthermore, for 
each value of riH shown, the lift coefficient for an open jet tunnel is 
lower than the infinite atmosphere value and increases monotonically with 
increasing ventilation coefficient until the lift coefficient exceeds the 
infinite atmosphere value. These calculation therefore indicate that for 
every finite value of tunnel depth to chord ratio, there exists a unique 
ventilation coefficient such that the lift coefficient equals the value 
of the lift coefficient in an infinite atmosphere. In the case of an 
infinite atmosphere, the lift coefficient is the same for all values of 
ventilation coefficient as is to be expected. 
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Table 12. Lift coefficient Cl^ vs. riH and C W for M = .85 and k = 0 


cw 

1 — 1 
II 

JC 

n H - 7.5 

o 

H 

II 

n H = ioo 

rifi — 1000 

0 

1.99486 (5) 1 

8.22740(5,10) 

8.98389(5) 

11.5788(5) 

11.8920(3) 

lcr 4 

1.99506(5) 

8.22744(5) 

8.98391(5) 

11.5788(4) 

11.8920(3) 

icr 2 

2.01449(5) 

8.23118(5) 

8.98630(5) 

11.5788(4) 

11.8920(3) 

i 

3.83187(5) 

8.57219(5) 

9.20734(5) 

11.5822(4) 

11.8920(3) 

10 4 

20.1757(5) 

11.8449(5) 

11.7403(5) 

11.7513(3) 

11.8952(3) 

h- 1 

o 

O' 

21.8952(5) 

12.2308(5) 

12.0980(5) 

11.9257(3) 

11.9243(3) 

loioo 

12.2351(5) 

12.0121(5) 

11.9292(3) 

11.9275(3) 

CO 


12 . 2442 (1) 2 




oo 


12 . 2351 (2-20) 2 





1 Numbers in parenthesis indicate the number of pressure basis functions used 
in the calculations. 

2 These results are those of Bland’s, which were programmed only for c^ = °°- 
3 Compare with nn = °°r CL a = 2tt /3 = 11.9275. 


Table 12 indicates that section coefficients are continuous with respect 
to ventilation coefficient for all values of ventilation coefficient. Since 
Ctyj can be any nonnegative real number, it is convenient for graphing purposes 
to map the domain of c w onto a finite interval. This can be accomplished 
by the transformation 

0^ = cot~ 1 c w (11-1) 

which is equivalent to restating the boundary condition (2-49) as 

p sin + Is cos e W = (11-2) 

dy 

77 

Then 0^ = 0 corresponds to a closed tunnel and 0 W = — corresponds to an 
open jet tunnel. Thus, we may call 0 W the ventilation angle . Graphs of 
section coefficients vs. ventilation angle are presented in Figures 14 to 16. 
Lift and moment are seen to be monotonically decreasing functions of venti- 
lation angle. The ventilation angle at which the lift coefficient equals 
the corresponding infinite atmosphere value differs from the ventilation 
angle at which the moment coefficient equals the corresponding infinite 
atmosphere value (alas) . The center of pressure, shown in Figure 16, is 
aft of the quarter chord for zero ventilation angle and is forward of the 
quarter chord for an open jet. 








Present method TWOBI 



Figure 16. Center of pressure vs. ventilation an^le for N 

k * 0 and n« “ 7*5 


M=i 




§12. Combined effects of depth to ch ord ratio an d wall ventilation 

This section presents predictions of the combined effects of tunnel 
depth to chord ratio and wall ventilation coefficient c w . Since the 
number of possible combinations of parameters can be quite large, we shall 
restrict the discussion to section coefficients defined for a limited domain 
of Mach number, frequency and mode shape; 

(1) M = 0, 1//2, 

(2) k = 0, .1, 

(3) vertical translation (h) , pitch rotation (a) , 

and to the doubly infinite domain of tunnel depth and wall porosity? 

( 4 ) 1 <_ tih < 00 r 

(5) 0 <_ c w < 00 . 

The resulting section coefficients can then be graphed as surfaces defined 
over a two dimensional region of tunnel depth and wall ventilation. It 
is convenient for graphing purposes to transform the supporting domain onto 
one of finite extent. This has already been partially accomplished in 
Section 11 by introducing the ventilation angle 0 W as defined by equation 
(11-1) . Similarly we define a tunnel depth angle 0 H by 

0H = cot -1 riH (12-1) 

7T 

which varies from 0 for infinite depth to — for zero depth. The depth 
angle is shown in Figure 17 as the angle between the line perpendicular 
to the airfoil passing through the midchord and the line from the midchord 
to the point of intersection with the tunnel wall of a line perpendicular 
to the airfoil passing through the trailing edge. 
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For present purposes, however, we shall consider only tunnels with depth 
greater than or equal to the airfoil chord. Thus, 

0 1 e H 1 (12-2) 

We can now represent the section coefficients for all possible values of 
ventilation coefficient and tunnel depth to chord ratio as surfaces with 
bounded support. 

CL a = CL a (9Hf 0 W ) » c Lh = c Lh (0 H> e w)f etc - (12-3) 

Figure 18 presents in compact form the complete range of values of 
lift coefficients Ci^ at M = 0 and k = 0 for all tunnel depths greater 
than or equal to the airfoil chord and for all possible values of wall 
ventilation coefficient. Several features of this CL a surface may be 
observed. The line 0 h = 0 corresponds to a tunnel of infinite depth and 
indicates a constant value of CL a = 2 tt for all values of ventilation co- 
efficient as it should, and is in keeping with the condition that the air- 
foil pressure not be affected by walls if they are infinitely far away. 

The line 0^=0 corresponds to a completely closed wall and indicates that 
Cl^ increases as the walls are brought closer together. The amount of 
increase in Cj^ with depth decreases with increasing ventilation so that 
beyond a certain value of 0^, CL a decreases as the walls are brought closer 
together. At the present time, we have not investigated if the trend 
reversal is reflected by experiment but this question is fundamental and 
bears on the validity of the boundary condition (2-49) . 

Figure 19 shows the surface for M = 1//2 and k = 0, and displays 

the same trends as Figure 18 for M = 0. Along the line 0 h = 0 (riH = 00 ) / 

2tt 

the value of Cl™ equals the expected value of — for all values of 0 W . 

This magnification with respect to Mach number is observed throughout the 

entire Ct surface, but the magnification factor is given by the well known 
1 

value — only for the infinite depth condition 0 h = 0. 

p 

Figures 20 and 21 show the surfaces for M = 0, k = 0 and M = l//2~, 

k = 0 respectively. Since the pitching moment is computed about the quarter 
chord, = 0 along the lines 0 H = 0 which correspond to the infinite 

depth case rijj - For closed tunnel walls, the moment coefficients increase 
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as the walls are brought closer together, the rate of increase being 
greater at the higher Mach number. For intermediate values of 0^ the 
increase with depth angle decreases until the moment coefficients decrease 
with increasing 0 H , similar to the behavior of the lift coefficient surfaces. 
Reversal occurs sooner at the higher Mach number. 

Figures 22 and 23 show the center of pressure surfaces for the two 
Mach numbers. In both figures the center of pressure is at the quarter 
chord for infinite depth (0 H = 0) , as is to be expected. For closed walls, 
the center of pressure moves aft as the walls are brought closer together, 
more so at the higher Mach number. Again, this trend reverses for venti- 
lated walls, more pronounced at the higher Mach number. 

Figures 24 to 27 show the surfaces representing magnitudes of lift 
and moment coefficients for vertical translation and pitch rotation about 
midchord at M = 0 and reduced frequency k = .1. Surfaces representing 
phase angles are shown in Figures 28 to 31. 



M =3 




M=i 



Figure 19. Lift coefficient C L v*. depth and ventilation for 
M = 1//? and k = 0 
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§13. _ Ef fpct^ pf yent ilati^n, on air foil- tunnel acoustic resonance 


Acoustic resonance between an airfoil and a ventilated wind tunnel 
will occur according to the frequency spectrum which is given by 

3^n 


kn ~ 


Mn H 


/ n " 1 , 2 , a . . 


(13-1) 


based on equation (2-53) . The values of reduced frequency at which reso- 
nance will occur depend upon Mach number/ tunnel depth to chord ratio and 
the ventilation coefficient. The effect of Mach number on the resonant 
frequencies consists in the factor 


£ 

M' 

which indicates that resonance can occur only for compressible flow and 
that the resonant frequencies approach zero as the Mach number approaches 
one. 


The effect of depth to chord ratio on the spectrum is most pronounced 
in the factor 


nH 

indicating that narrow tunnels have higher resonant frequencies whereas 
deeper tunnels have lower resonant frequencies/ as is to be expected. 
Reversing the argument/ it may be noted that for a given Mach number and 
reduced frequency, there exist infinitely many values of depth to chord 
ratio at which acoustic resonance will occur: 

$A n 

nH(n) = 1ST' n = 1 ' 2 '"* 


(13-2) 


The effect of ventilation coefficient cw on the frequency spectrum is 
combined with the depth to chord ratio rin through the relation 

Cm 

An = * n (— >• (13-3) 

•H 

Referring to Figure 7 , the effect of ventilation is to shift the spectrum 
within the bounds 


closed - tt (n— “) A n <_ im ~ open 

so that increasing the amount of ventilation increases the resonant frequencies. 
The maximum fractional change in frequency that can be brought about by 
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ventilation alone is to double the fundamental frequency in going from 
a closed wall to an open jet- This is illustrated by Table 13. 

Table 13. Resonant frequencies vs. ventilation 

y/J 

coefficient for M = and Tin = 10 



GW = 00 
(closed) 

c w ~ 1 
(ventilated) 

C W = 0 

(open) 

k l 

.090690 

.165282 

.181380 

k 2 

.181380 

.332586 

.362760 

k 3 

.272070 

.460648 

.544140 


The effect of acoustic resonance is to cause the pressure to vanish 

at the resonant frequency. This is illustrated by Figure 32 for a flat 

/3 

plate oscillating about the midchord at M = — , with three values of venti- 
lation coefficient corresponding to a closed wall (c^ = °°) , a ventilated 
wall (c w = 1) , and an open jet (c w = 0) . For comparison, the value of 
|C L J at M = 0 is shown also since resonance cannot occur in incompressible 
flow. The comparison is striking. Whereas the behavior of vs. reduced 

frequency is smooth for M = 0 with the three curves for c^ = Of If and 00 

• . 1 1 1 / 3 " * 

merging as the frequency increases, the behavior of | | for M = is 

quite different. The closed wall condition begins with a relatively large 

value of CL a at k = 0 and drops to 0 at resonance very abruptly, increasing 

for values of reduced frequency beyond resonance to a maximum value of 

approximately 60% of its zero frequency value, then dropping to 0 again 

at the second resonant frequency, and so on. Similar behavior is evidenced 

by the ventilated wall conditions, beginning with a lower value of C L a at 

k = 0 and displaying higher resonant frequencies. In all cases the drop to 

CL a = 0 at resonance is more abrupt below resonance than above, and the 

abruptness appears to decrease slightly at the higher resonant frequencies. 

The phase angle is shown in Figure 33 with and without resonance. No 

resonance occurs for incompressible flow. In all cases a 90 degree phase 

lag is predicted at resonance, dropping very abruptly above resonance until 

the next resonant frequency, at which point a 90 degree phase lag occurs 

again, and so on. 
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Figure 32. Effect of ventilation on I C T 

a 

resonance effects 


vs. k with acoustic 
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§ 14. Co ncluding remarks 

The numerical calculation of unsteady airloads on thin airfoils in subsonic 
ventilated wind tunnels has been accomplished using Bland's kernel (2-50). 
Additionally, we have rigorously proved that the collocation method of 
solving Bland's integral equation (2-57) converges to the mathematically 
exact solution, and have established a new three-way equivalence between 
collocation, least squares and Galerkin's method whenever the collocation 
points are chosen as the nodes of the quadrature rule used for Galerkin’s 
method. This convergence behavior has been demonstrated with the TWODI 
program. We point out that the convergence proof, given for the first time 
by this work, applies to an arbitrary kernel whose dominant singularity 
is of Cauchy type, and is thereby of a general nature. Furthermore the 
method of proof, based on converting the integral equation to one of the 
second kind, opens up the way to methods of solution and error estimates 
not otherwise available for equations of the first kind. 

Results from the computer program have been compared with exact closed 
form solutions in special cases. Using NP <_ 10, six decimal accuracy is 
obtained for steady flow and three decimal accuracy or better is attained 
for unsteady flow. New results are presented showing the effect of wall 
ventilation and depth to chord ratio on section coefficients, and on acoustic 
resonance between the airfoil and the tunnel walls. While these results 
should improve the confidence and precision of wind tunnel testing, it 
would be desirable to compare the predictions of the TWODI program with 
known experimental data, and with future experimental data for unsteady 
flow in ventilated tunnels as they become available. 

Although the TWODI program is already sufficiently accurate for most 
engineering purposes, large amounts of computer time are required with the 
pilot version for very deep tunnels, particularly in unsteady flow. This 
time is largely expended in subroutine SUM, and more efficient numerical 
algorithms would probably alleviate this difficulty. At higher frequencies, 
we have observed a deterioration of the rate of solution convergence with 
respect to the number of pressure basis functions. Preliminary indications 



are that the kernel is evaluated accurately but may not be integrated 
accurately for high frequencies. Since Hsu's interdigitation procedure 
restricts the number of quadrature points to equal the number of pressure 
basis functions , it is obvious that for a fixed value of NP, there is a 
frequency k at which the factor 

e -ikx 

in the continuous part of the kernel will render the Jacobi-Gaussian 
quadrature inaccurate. (This provides additional incentive for improving 
the efficiency of computing the infinite series in subroutine SUM.) While 
the physical justification of linearization fails at sufficiently high 
frequency [89,1948], a high frequency capability enables the complete 
frequency spectra to be calculated, from which point the response to 
arbitrary time dependent excitation can follow. It is therefore potentially 
valuable to study the problem of numerically computing 

( x_ £) f (x-£) d£ 

accurately for large k and continuous functions f, using a small number of 
function evaluations. 

Although the convergence proof applies to discontinuous downwash 
functions as well, the rate of convergence for airfoils with flaps may 
be expected to be weaker owing to a form of Gibb's phenomenon for general- 
ized Fourier series. Therefore, as a practical matter, it would be 
desirable to extend the solution method to permit the efficient analysis 
of multi-control airfoil configurations. 

Bland's kernel is based upon a boundary condition which is phenomeno- 
logically approximate whenever the tunnel walls are other than fully closed 
or an open jet. Removal of this restriction awaits a more rational theory 
for ventilated wind tunnel walls. 
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COMPARISON MITH THE KUESSNER-SCHKSR2 SOLUTION USING 
77/65/14. 39.15.51. ?fiGE 
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* CSLCULftTION OF THO-D 1NENS I ONAL- UNSTEADY SIRLOS3S * 
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* SUBSONIC COMPRESSIBLE FLOH WIND TUNNEL * 
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CASE 


4 
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HACH NO. FREQUENCY 

0.000000 1.000000 


ETR-SUB-K 

380.000 


c-sub-h 
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CORPflRlSOS iilTH THE KUESSNER-SCHH8RZ SOLUTION U8IN8 ETRH = 388 


77/85/14. 




CASE 1 OF 1 
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FREQUENCY ETfl-SUB-H 
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COMPARISON WITH TH£ KUESSNES-SCHSiRRZ SOLUTION USING ETflN 
77/06/14. 10.35.13. PAGE 3 


DOHNNASH NODE 4 
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COHPhR ISO fi WITH THE KliESS^ER-SCHHRRZ SOLUT i ON USING ETnH = 300 
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4.48792 

-47.i79 

i 

-.380000 

■ -J V V V I s 0 

3 « it*o**t) 


-2.33527 

4.39160 

-41.360 

i 

-.200900 

< 4 V V U V V 

3.42353 


- £ * & 3 9 0 $> 

4.32264 

-37.627 

i 

-.100009 

.458000 

0 K ■ J li d 

yiyyi’r't 


-2.38253 

4.26092 

-34.007 

i 

8.000000 

.588800 

o s o j d o 


-2.15564 

4.19841 

-30.959 

i 

.108880 

.550808 

-J ■ 'J * i tJ u 


-1.95848 

4.18453 

-28.372 

i 

.208080 

.680880 

3.58626 


-1.76132 

3.99544 

“tfbi i57 

i 

.380000 

.650880 

3.51642 


- i i 53352 

3.85656 

-24.244 

i 

.400000 

.709088 

3.33913 


-1.41348 

3.63130 

-22.579 

i 

.500000 

.750800 

3.22903 


-1.24716 

3.46151 

-21.118 

i 

.600000 

Q fi ft ft ft ft 
> 0 y v v y v 

2.39694 


-1.08055 

3.18579 

-19.827 

i 

.700080 

.350008 

2.63663 


-.308201 

2.83598 

-18.678 

1 

.380000 

.900080 

2.26536 


-.720862 

2.37776 

-17.643 

4 

* 

.980000 

.950008 

1.65186 


-.496299 

1.72473 

-16.720 

i 

1.000000 

i. 808000 

8. 


0. 

0. 

9.000 

2 

-.900000 

.050000 

. 1166 


13.9362 

26.1412 

i47. 784 

2 

-.380000 

.190080 

-17.2555 


4.78367 

17.3079 

164.48? 

2 

-.780000 

.158008 

-15.1143 


-.675856E-92 

15.1143 

-179.974 

2 

-.600000 

.299800 

-13.7333 


-3.20222 

14.1917 

-166.375 

2 

-.500000 

.250908 

-12.6514 


-5.53986 

13.8111 

-156.352 

2 

-.480000 

. v 9 0 8 3 0 

-11.7045 


-7. 32451 

13.8074 

-147.962 

2 

-.300000 

.358089 

-10.8222 


-8.70957 

13.8916 

-141.173 

2 

-.200000 

.488090 

-9.97075 


-3.78131 

13.3679 

-135.548 

2 

-.100000 

.459800 

-3.13379 


-10.5941 

13.3379 

-130.767 

2 

8.000808 

.508088 

-3.30354 


-11.1789 

13.9254 

-126.605 

2 

.180000 

.550000 

-7.47705 


-11.5560 

13.7640 

-122.904 

2 

.280000 

.690000 

-6.65414 


-11.7358 

13.4909 

-119.553 

2 

.300000 

.650000 

-5.33641 


-11.7203 

13.0931 

-116.472 
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COMPARISON HITH THE KUESSNER-SCHHARZ SOLUTION USING ETftH = 300 
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.400000 .700000 

.500000 .750000 
.600000 .800000 
.700000 .850080 
.800000 .988000 
.900000 .950000 
1.000000 1.000000 
-.980000 .050000 
-.800000 .100000 
-.700000 .150000 
-.600000 .200000 
-.500000 .250000 
-.400000 .300000 
-.300000 .350000 

-.200000 .400000 
-.100000 .450000 
0.000000 .500000 

.100000 .550000 
.200000 .600000 
.300000 .650000 
.400000 .700000 
.500000 .750000 
.600000 .300000 
.700000 .850000 
.800000 .900000 
.900000 .950830 
1.080000 1.800000 
-.900000 .053000 
4 -.800000 .130000 

4 -.700000 .150000 

4 -.600000 .200000 

4 -.500000 .250000 

4 -.400000 .300000 

4 -.300000 .350000 

4 -.200000 .400300 

4 -.100000 .450000 


-5.02650 

-11.5038 

-4.22769 

-11.0703 

-3.44334 

-10.3890 

-2.67593 

-9.40307 

-1.92400 

-7.99723 

-1.16821 

-5.87344 

0 . 

0 . 

39.8107 

5.58926 

13.9831 

i5 . 3537 

7.84448 

19.7334 

.137027 

21.7247 

-5.75172 

22.1258 

-10.4495 

21.4587 

-14.2549 

20.0112 

-17.3276 

17.9838 

-18.7591 

15.5292 

-21.6039 

12.7733 

-22.8940 

9.82604 

-23.6464 

6.73949 

-23.8664 

3.76377 

-23.5477 

.353010 

-22.6696 

-1.82716 

-21.1890 

-4.13739 

-19.0199 

-5.89500 

-15.9758 

-6.81261 

-11.5389 

-6.30739 

0 . 

0 . 

-57.9898 

-10.7160 

-10.3804 

-23.0932 

11.3234 

-25.0051 

22.4443 

-22.5341 

27.4200 

-17.9662 

23.2543 

-12.5350 

26.1017 

-7.01813 

21.7403 

-1.94244 

15.763? 

2.32700 


12.5540 

-113.603 

11.8501 

-110.902 

10.9448 

-108.337 

9.77642 

-105.885 

8.22541 

-103.527 

'5.98849 

-101.249 

0 . 

0.000 

40.2011 

7.992 

24.4150 

38.966 

21.2819 

68.371 

21.7251 

89.639 

22.8612 

104.572 

23.8677 

115.964 

24.5693 

125.464 

24.9732 

133.935 

25.1312 

141.835 

25.0975 

149.406 

24.9136 

156.771 

24.6018 

163.980 

24.1614 

171.038 

23.5632 

177.925 

22.7431 

-175.392 

21.5893 

-168.950 

19.9125 

-162.780 

17.3677 

-156.905 

13.1503 

-i51. 338 

0 . 

0.000 

58.9716 

-169.530 

25.3189 

-114.204 

27.4495 

-65.637 

31.8046 

-45.114 

32.781? 

-33.234 

30.9105 

-23.924 

27.028? 

-15.050 

21.3269 

-5.106 

15.9346 

8.397 
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COMPARISON HITH THE KUESSNSR-SCHNARZ 
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SOLUTION USING ETAH =300 
PAGE 6 


4 

0.000000 

• d 6960U 

3.66429 

5.5420? 

10.2852 

32.605 

4 

.100000 

.550080 

.879524 

7.5433? 

7.59943 

33.354 

4 

.200000 

.609000 

-7.13029 

3.27314 

10.9545 

130.955 

4 

.300000 

.650890 

-15.1117 

7.72111 

16.9699 

152.936 

4 

.400000 

,700080 

-22.4933 

5.97701 

23.2792 

165.123 

4 

.500000 

.750000 

-23.3890 

3.21639 

29.0674 

173.648 

4 

.600000 

.300000 

-33.7529 

-.272434 

33.7540 

-179.538 

4 

.700000 

.350000 

-36.4031 

-4.03566 

36.6261 

-173.674 

4 

•806089 

.900800 

-35.7343 

-7.33642 

36,5286 

-168.414 

4 

.900000 

.950900 

-29.6948 

-8.75209 

30.9577 

-163.578 

4 

1.000000 1.000000 

0. 

0. 

0. 

0.000 

5 

-.900000 

. 656636 

31.1405 

40.2407 

50.882? 

52.265 

5 

-.800000 

.108000 

-44,3302 

41.2824 

61.0528 

137.455 

5 

-.700000 

.159006 

-67.2518 

3i.9825 

74.4693 

154.566 

5 

-.600800“ 

'.288tf0‘8 

-fa?. 8753 

2873623 

7 0 - 3 6 i 3 

163. 301 

5 

-.500000 

.250000 

-57.9520 

9.71526 

58.7607 

170.483 

5 

-.400000 

.308008 

-*>3, 2ii-i4 

1.55265 

43.2312 

177.942 

5 

-.300000 

.350008 

-2? .1849 

—3. 55^ i 4 

2? . 416 0 

-172.554 

5 

-.280000 

.400908 

-12.2385 

-5.60663 

13.5071 

-155.475 

5 

-.100000 

.458000 

-.166123 

-4.99352 

4.99628 

-91.905 

5 

0.000000 

.508080 

8.07256 

-2.34813 

8.40714 

-16.219 

5 

.109000 

.550809 

11.7605 

1.53832 

if. 8607 

7.455 

5 

.200090 

.606009 

10.6249 

5.79789 

12.1039 

28.621 

5 

.390000. 

.650008 

A 77“>C7 

7 « l » £. C* 1 

9.55125 

10. 6773' 

63.449 

5 

.400009 

.700000 

-5.29617 

11.9839 

13.1021 

113.843 

5 

.500000 

.750800 

-18.6395 ' 

12.4243 

22.4008 

146.314 

5 

.600000 

.800000 

— -5 3 . { 6 4 ( 

10.4476 

35.3441 

162.807 

5 

.700000 

.350009 

-43.3910 

6.02970 

48.7652 

i?2.89? 

5 

.800000 

.298000 

-58.8496 

-.168415 

58.8499 

-173.836 

5 

.900009 

.950008 

-57.9353 

-6. 05846 

53.3010 

-174.035 

5 

i . oodooo 

1.000008 

0. 

0. 

0. 

0.000 
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SECTION COEFFICIENTS 


HOKE 

COEFF 

REAL 

IMAGINARY 

MAGNITUIIE 

PHASE ANI 

1 

LIFT 

2.50990 

-3.38758 

4.21608 

-53.465 

1 

MOMENT 

1.56898 

-.519905E-03 

1.56898 

019 

1 

X-CP 

.471642 

.298341 

.558401 

32.368 

2 

LIFT 

-9.91313 

-5.01876 

11.1112 

-153.148 

2 

HQHENT 

-.786022 

-6.27628 

6.32531 

-97.138 

2 

X-CP 

.568254 

.472004 

.738715 

39. 714 

3 

LIFT 

-6.76970 

7.52965 

10.1254 

131.358 

3 

MOMENT 

-13.3370 

.224304E-02 

13.3370 

i79.990 

3 

X-C? 

i. 13081 

.979352 

1.49595 

40.895 

» 

LIFT 

-13.5648 

-3.77175 

14.0794 

-164.461 

4 

MOMENT 

-6.27720 

-.119256E-01 

6.27721 

-179.891 

4 

X-CP 

.679774 

-.118621 

.690046 

-3.839 

5 

LIFT 

-13.5511 

3.80183 

16.1588 

146.995 

5 

MOMENT 

-18.8494 

.920512E-03 

18.8434 

179.997 

5 

X-CP 

1.22830 

.635361 

i. 38239 

27.351 
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ROW i 
ROR 2 
RON 3 
ROH 4 
RU9 5 


COflPfiRiSQfi aiTH THE KUES3NES-SCHURR2 SOLUTION USING ETftH = 300 
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GENERALIZES 

fl'ERODYUUHIC FORCES 


2.50990 

-3.38758 

-9.91313 

-5.01876 

-6.76970 

7.52965 

“13. 5648 

-3.77175 

-13.5511 

8.80183 



-1.88183 

6.774i2 

18.2542 

-2.51504 

-13.1345 

-15. 0548 

14.5752 

7.51964 

-10.5964 

-17.6018 



-1.25484 

U! 1 1 UJX 

-15.1216 

15.0644 

39.2777 

2.51067 

-40.7266 

-7.48823 

35.7409 

17.5809 



1.25551 

6.77321 

13.5613 

-15.0695 

-37.7213 

A A* iV'‘kv 

iili i/tw 

79.1556 

—5 . U5is 0 

-74.1779 

-17.5537 



-1.25751 

-6.77809 

-13.5587 

15.0776 

36.6959 

-10. 0346" 

-78.1295 

i?.5?89 

124.989 

5.02429 
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